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ABSTRACT 

A viable alternative to the dark energy as a solution of the cosmic speed up 
problem is represented by Extended Theories of Gravity. Should this be indeed the 
case, there will be an impact not only on cosmological scales, but also at any scale, 
from the Solar System to extragalactic ones. In particular, the gravitational potential 
can be different from the Newtonian one commonly adopted when computing the 
circular velocity fitted to spiral galaxies rotation curves. Phenomenologically modelling 
the modified point mass potential as the sum of a Newtonian and a Yukawa- like 
correction, we simulate observed rotation curves for a spiral galaxy described as the 
sum of an exponential disc and a Navarro - Frenk - White (NFW) dark matter halo. We 
then fit these curves assuming parameterized halo models (cither with an inner cusp or 
a core) and using the Newtonian potential to estimate the theoretical rotation curve. 
Such a study allows us to investigate the bias on the disc and halo model parameters 
induced by the systematic error induced by forcing the gravity theory to be Newtonian 
when it is not. As a general result, we find that both the halo scale length and virial 
mass are significantly overestimated, while the dark matter mass fraction within the 
disc optical radius is typically underestimated. Moreover, should the Yukawa scale 
length be smaller than the disc half mass radius, then the logarithmic slope of the 
halo density profile would turn out to be shallower than the NFW one. Finally, cored 
models are able to fit quite well the simulated rotation curves, provided the disc 
mass is biased high in agreement with the results in literature, favoring cored haloes 
and maximal discs. Such results make us argue that the cusp/core controversy could 
actually be the outcome of an incorrect assumption about which theory of gravity 
must actually be used in computing the theoretical circular velocity. 
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1 INTRODUCTION 

The picture of a spatially fiat universe undergoing accel- 
erated expansion is the nowadays accepted view of our 
cosmo. According to the s u ccessful concordance ACDM 
model (|Carroll et al. I 1 19921 ; ISahni fc Starobinski I |2000| ). 
there are two main ingredients in this scenario, namely 
dark matter (accounting for the clustering of the struc- 
tures we observe) and the cosmological constant A (dom- 
inating the energy budget and driving the cosmic speed 
up). The anisotropy spectrum of c osmic microwave back - 
ground (|de Bernardis et al. 1 12000| ; iKomatsu et al. I l20TOl h 
the galaxy power spectrum with the imprinted baryon 



acous tic oscillations |Eisenstein et al.~ll2005l ; IPercival et al. I 
|2010T) and the Hubb l e diagram of Typ e la Supernovae 
(IKowalski et al. I l2008l ; iHicken et al. II2OO9] ) represent an in- 
complete list of the wide amount of data this model is able 
to excellently reproduce in a single scenario. 

From a theoretical point of view, however, the ACDM 
is far to be satisfactory. First, the A term is ~ 120 orders 
of magnitude larger than what expected from quantum field 
theory, while the ratio between its energy density and the 
matter one is coincidentally of order unity just today while 
it should have been much smaller or much larger than 1 
over the rest of the universe history. Motivated by these 
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unpleasing shortcomings, a plethora of models giving rise 
to a varying A - like term has been proposed mainly based 
on a scalar field evolving under the influence of its own self 
interaction potential. Needless to say, the absence of any 
candidate for this scalar field and the full arbitrariness in 
the choice of the potential are serious drawbacks of these 
models which therefore represent only a way to change the 
problems without actually solving them. 

On galactic scales, the A term gives a negligible con- 
tribution to the gravitational potential so that the classical 
Newtonian theory is usually adopted. While the evolution of 
the universe is driven by the cosmological constant, the for- 
mation of structure is mainly determined by the dark matter 
(DM) which provide the potential wells where baryons col- 
lapse to originate the visible component of galaxies. Numer- 
ical simulations allow to follow this process predicting the 
structure of DM haloes. Surprisingly, the theoretical expec- 
tations are not in agreement with observations on galactic 
scales. In particular, the density profile of DM haloes is ex- 
pected to follow a double power - law with p tx x~ a (l+x) 3 ~ a 
with x — r/r s , r s a characteristic length scale and a giving 
the logarithmic slope in the inner regions. Although the de- 
bate on what the precise value of a is remains open, what is 
granted is that a should definitel y be positive with a — 1 for 
the most popular NFW model ^Navarro et al. 1 1 19971 ). The 
DM haloes are thus described by cusped models or, in other 
words, the logarithmic slope 7 = din p/d\xir never vanishes. 
On the contrary, rotation curves of low surface brightness 
galaxies (LSB) are definitely better fitted by cored models, 
i.e. 7 = in the inner regi ons, such as the pseudo - iso thermal 
sphere, p oc 1/(1 + x 2 ) (jBinnev fc Tremahie lll987n or the 
Burkert model, p oc (1 + z)"^! + x 2 )' 1 (|Burkertl Il995l ; 



ISalucci fc Burkert llioooh . As well reviewed in de Blok et al 
(2010), cored models turn out to be statistically preferred 
over cusped ones from the fit to large samples of LSB rota- 
tion curves. Moreover, for those galaxies where a statistically 
acceptable fit for a cusp model is obtained, it turns out that 
the concentration parameter (defined later) is significantly 
larger than wh at expected from numerical simulations (see 
Ide Blokll2010l and refs. therein for further details). 

From the above picture, it is clear that, although ob- 
servationally successful on cosmological scales, the ACDM 
model is far to be free of problems, especially if one also 
remember that there is up to now no laboratory final ev- 
idence for any of the many particles candidate to the role 
of DM. It is therefore worth going beyond the usual view 
and look at a radical revision of the underlying scheme. To 
this end, one can consider cosmic speed up as the first evi- 
dence of a breakdown of General Relativity (GR) as we know 
it. Rather than being due to a new actor on the scene, the 
accelerated expansion can be a consequence of gravity work- 
ing in a different way than GR predicts. This consideration 
has attracted most interest towards braneworl d - like mod- 
els such as the five dimension s DGP models (jDvali et al. I 
l200d : iLue fc Starkman I I2003T ) or fourth order theories of 
gravity, where the GR Einstein - Hilbert Lagrangian is gen- 
eralized by the introduction of a function f(R) of the 
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mological scales, one has to check whether the gravitational 
potential on galactic scale is. Most of modified theories in- 
deed induce negligible changes to the gravitational potential 
on Sola r System sc ale in order to pass the classical tests of 
gravity (|Will Ill993t) . but this does not prevent to have signif- 
icant deviations from the Newtonian potential on the much 
larger scale of galaxies where no direct experimental test is 
available. 

Although modified gravity theories are investigated at 
cosmological scales as alternatives to dark energy, we are 
here interested in whether they can also impact the estimate 
of dark matter properties on galactic scales. Should indeed 
the gravitational potential be modified, the computation of 
the rotation curve must be done in this modified framework. 
On the contrary, one typically assumes that Newtonian me- 
chanics holds and then constrains the halo model parameters 
by fitting the theoretical rotation curve to the observed one. 
We are here interested in investigating whether this incor- 
rect procedure biases in a significant way the determination 
of the halo parameters and whether such a bias can explain 
the inconsistencies among theoretical expectations and ob- 
servations. In a sense, we are wondering whether modified 
gravity could be a possible way to solve the cusp/core and 
similar problems of the DM scenario. 

The plan of the paper is as follows. In Section 2, we 
present the modified gravitational potential used and de- 
tail the derivation of the rotation curve. It is worth noticing 
that Yukawa-like corrections emerges in any analytical f(R)- 
gravity model, except f(R) = R, where R is the Ricci scalar 
adopted in the Hilbert-Einstein action. Section 3 describes 
how we estimate the bias induced on the halo model pa- 
rameters by fitting data with the Newtonian potential while 
the underlying theory of gravity is non-Newtonian. The re- 
sults of this analysis are discussed in Sections 4 and 5, while 
Section 6 summarizes our conclusions. 



2 YUKAWA-LIKE GRAVITATIONAL 
POTENTIALS 

As it is well known, Newtonian mechanics is the low energy 
limit of GR. Indeed, looking for a stationary and spheri- 
cally symmetric solution of the Einstein equations gives the 
Schwarzschild metric whose tt component gives to the New- 
tonian 1/r gravitational potential in the weak field limit. 
Modifying GR leads to modified field equations hence to a 
different solution and potential in the low energy limit. As 
such, we should first choose a modified theory of gravity to 
finally get the modified potential. 

An interesting example for its application to cos- 
mology is provided by f(R) theories of gravity. Writ- 
ing the metric in the weak field limit as 

ds 2 = -[1 - 2A(r)]dt 2 + [1 + 2B(r)]dr 2 + r 2 dQ. 2 , 

and using the conformal equivalence with scalar - 
field theories, Faulkner et al. (2007) have shown 
that : 



A(r) = A(f) + 



(1) 



It is worth stressing that, notwithstanding which is the cor- 
rect modified gravity theory, should GR be incorrect on cos- 



where tilted quantities are evaluated in the Einstein 
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frame (with f = x and \ = f'(R)) and ip(r) is the 
scalar field coupled to matter determined by : 



l__d ( 2 chp\ _ dV ef f{1>,f) 



f 2 df 



df J 



(2) 



The effective potential V e // is determined from both 
the functional expression for f(R) and the local mass 
density p(r) as : 



V eff = VW + X ~ 1/2 p(r) 



(3) 



with V(V0 given by Eq.(6) in Faulkner et al. (2007) 
and p = X~ S ^ 2 P- For the particular case of a uniform 
sphere of mass M c and radius R c and a quadratic 
potential, one finally gets : 

(f>(r) oc - [1 + (1/3) exp (-m^pr)] 

for the gravitational potential outside R c . For the 
general case, depending on the choice of f{R), the 
chamaleon effect (see, e.g., Faulkner et al. 2007 and 
refs. therein) can take place leading to the same 
above potential but with a different factor A/3 in- 
stead of (1/3), with A << 1 depending on both the 
source mass and the local density which is embedded 
in. 

Motivated by this consideration, we will postulate that 
the gravitational potential generated by a point mass m is : 



(r) 



Gm 
'(1 + S)r 



1 + 5 exp ( — — 



(4) 



where (5, A) depend on the parameters of the theory. It 
is worth noticing that a Yukawa- like correction has been 
invoked several times in the past. For instance, Sanders 
(1984) showed that the observed fiat rotation curves of spi- 
ral galaxies may be well fitted by this model with no need 
for dark haloes provided 5 < and A is adjusted on a case - 
by -case basis. As discussed above, Yukawa- like corrections 
have been obtained, as a ge neral feature, in the frame work 
of /(_R)-theories of gravitjo (jCapozziello et al~ll2009ah and 
successfully applied to clu sters of galaxies setting 5 = 1/3 
jCapozziello et al. Il2009bh . In general, one can relate the 
length scale A to the mass of the effective scalar field intro- 
duced by the Extended Theory of Gravitjjf]. The larger is 
the mass, the smaller will be A and the faster will be the 
exponential decay of the correction, i.e., the larger is the 



1 As hinted at above, for f(R) theories showing the 
chamaleon effect, 5 should be a function of the source 
mass m. We will therefore implicitly assume that all the 
stars in a galaxy have the same mass. Needless to say, 
this is far to be true, but, as far as the value of 5 does not 
change too much with m, the impact of this simplifying 
assumption can be neglected. Alternatively, one should 
introduce an average over the stellar mass function, but 
this latter quantity is largely unknown so that we prefer 
not to make any arbitrary assumption on its functional 
form and fully neglect the dependence of <5 on m. 

2 Referring to /(i?)-gravity, we prefer to deal with "Extended 
Theories of Gravity" and not modified theories of gravity since 
these theories are nothing else but a straightforward extension of 
GR w here the considered action is a generic function of the Ricci 
scala r (Capozzicl lo fc Francaviglia 1 I200& ICapozziello fc Faraoni I 
1201(f) . 



mass, the quicker is the recovering of the classical dynam- 
icfl Eq.Q then gives us the opportunity to investigate in 
a simple and unified way the impact of a large class of mod- 
ified gravity theories, among these the Extended Theories, 
since other details do not have any impact on the galactic 
scales we are interested in. 

Eq.Q is our starting point for the computation of 
the rotation curve of an extended system. To this end, we 
first remember that, in the Newtonian gravity framework, 
the circular velocity in the equatorial plane is given by 
Vc(R) — Rd&/dR\, = o, with $ the total gravitational poten- 
tial. Thanks to the superposition principle and the linearity 
of the point mass potential on the mass m, this latter is com- 
puted by adding the contribution from infinitesimally small 
mass elements and then transforming the sum into an inte- 
gral over the mass distribution. For a spherically symmetric 
body, one can simplify this procedure invoking the Gauss 
theorem to find out that the usual result v c (r) = GM(r)/r 
with M(r) the total mass within r. However, because of the 
Yukawa -like correction, the Gauss theorem does not apply 
anymore and hence we must generalize the derivation of the 
gravitational potential- Alternatively, one can also remem- 
ber that v c (R) = RF{R,z = 0) being F the total gravita- 
tional force and R a radial coordinate. This is the starting 
point adopted in Cardone et al. (2010) where a general ex- 
pression has been derived for the case of a generic potential 
giving rise to a separable force, i.e. : 



F P (p,r) 



—2-Uwfr{V) 



(5) 



with fi = itl/Mq, rj = r/r s and (Me,r s ) the Solar mass 
and a characteristic length of the problem. In our case, it is 
— 1 and : 



1 + 



exp (~r]/rjx) 
(l + 5)?7 2 



(6) 



with rjx — X/r s . Using cylindrical coordinates (R,9,z) and 
the corresponding dimensionless variables (77, 8, £) (with £ = 
z/r s ), the total force then reads : 



Fir) 



GpQTs 



n'dr,' / d(' / f r (A)p( v ',e',(')de' (7) 

'0 J-oo JO 



with p = p/po, po a reference density, and we have defined 



A = [n + V 2 - 2W cos {6 - 6') + (C - C') 2 ] 



,x2-|l/2 



(8) 



Since we will be interested in axisymmetric systems, we can 
set p — p{r],C,). Moreover, the systems of interest here are 
spiral galaxies which we will model as the sum of an infinites- 
imally thin disc and a spherical halo, so that a convenient 
choice for the scaling radius r s will be the disc scale length 



3 Note that the factor 1/(1 + 5) has been explicitly introduced in 
order to recover the correct Newtonian potential for A — > oo. 

4 The non-validity of the Gauss theorem is not a shortcoming 
since, as discussed in Capozziello et al. (2007), physical conserva- 
tion laws are guaranteed by the Bianchi identities that must hold 
in any modified theories of gravity. 
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Rd- Under these assumptions, the rotation curve may then 
be evaluated as : 



iii. Fit the simulated curve with the same disc + halo 
model, but using the Newtonian theory to compute v c (R). 



GpoRlv 



vi(R) = 



1 + 5 



r,'dr,' / p(r/,C'R' / MA o )d0' (9) 

>0 J-oo JO 



with 



A O = A(0 = C = O) = [r, 2 + V ' 2 -2 V r,'cose' + C' 2 ] 1/2 .(10) 

Inserting Eq.© into Eq.([9| gives rise to an integral which 
has to be evaluated numerically even for the spherically sym- 
metric case. It is, however, clear that the total rotation curve 
may be splitted in the sum of the usual Newtonian one and 
a corrective term disappearing for A — > oo, i.e. when the ex- 
tended gravity has no deviations from GR on galactic scales. 



3 ESTIMATING THE BIAS 

Let us assume that a spiral galaxy can be modelled as the 
sum of an infinitesimally thin disc and a spherical halo and 
denote with p the halo model parameters. We can then write 
the total rotation curve as : 



v 2 (R, M d , Pi 



= v} N {R,M d ) +vl N (R,pi) 
+ v 2 dY (R,M d ) + v 2 hY (R, Pl ) 



where M d is the disc mass, the labels d and h denote disc and 
halo related quantities, while N and Y refer to the Newto- 
nian and Yukawa - like contributions. These latter terms fade 
away for r >> A so that the outer rotation curve is likely the 
same as the Newtonian one. On the other hand, in the in- 
ner region, the two curves may differ more or less depending 
on the value of X/Rd- It is worth wondering whether such 
a difference may be compensated by adjusting the model 
parameters. That is to say, we are looking for a new set 
(M' d ,p') such that: 



2 



(R,Md,p) = v 2 dN (R,M d ) + v 2 hN (R,p') . 



Formally, this problem could be solved explicitly writing 
down the above relation for M + 1 values of r with M 
the number of halo parameters and then checking that the 
matching between the two curves is reasonably good (if not 
exact) along the full radial range. Actually, such a procedure 
is far from being ideal since it introduces a dependence of 
the results on the values of r chosen. Moreover, we do not 
need to exactly match the two curves, but only find (M' d , p') 
in such a way that the two curves trace each other within 
the typical observational uncertainties. 

In order to find (M' d ,p') taking care of typical observa- 
tions, we therefore adopt the procedure sketched below : 

i. Compute the theoretical circular velocity for input 
model parameters (Md,p) using the exact expression. 

ii. Generate a simulated rotation curve by sampling the 
above v c and adding noise to the extracted points. 



iv. Compare the output parameters (M' d ,p') with the 
input ones (M d ,p) as function of the normalized scale 
length r)\ — X/Rd of the modified gravity theory and other 
quantities of interest. 



In the following, we describe in more details steps i. and ii., 
while Sections 4 and 5 are devoted to discussing the results 
from a large sample of simulated curves. 



3.1 Modeling spiral galaxies 

As yet said above, we will model a spiral galaxy as the sum 
of a thick disc and a spherical halo. For the disc, we adopt 
a double exponential disc so that the density profile reads : 



Pa(R, 



M d 



■ exp 



R_ 

Rd 



(11) 



where (Md,Rd,Zd) are the disc mass, lengthscale and 
heightscale. The Newtonian rotation curve cannot be com- 
puted analytically, but, if Zd « Rd (as we will assume), it 
is well ap proximated by t he formula for the infinite simally 
thin disc ^Freeman Ill97d : iBinnev fc Tremaine Ifl987l ) : 

v 2 dN {R) = {2GM d /Rd)y 2 [I (y)K (y) - h(y)Ki(y)] (12) 

with y — R/2Rd and I n (y), K n (y) the modified Bessel func- 
tions of order n of the first and second kind, respectively. 

While the observed photometry motivates the use of 
the exponential profile for the disc, the choice of the dark 
halo model is not trivial. Numerical simulations of struc- 
ture formation are typically invoked as a direct evidence 
favouring the use of the NFW density law or its variants. 
However, the NFW model is the outcome of DM only sim- 
ulations performed in a Newtonian framework, while here 
we are working in a modified gravity theory. In principle, 
one should therefore rely on the results of simulations which 
include both the effect of the different potential and the im- 
pact on the evolution of structure due to deviations from 
GR. To this end, one has first to specify which is the mod- 
ified gravity theory one is considering, i.e., explicitly write 
down the gravi t y Lag rangian. In the case of /(-R)-gravity, 
ISchmidt et al. I ( 2009h have shown that, provided f(R) sat- 
isfies the constraints from the Solar System, the halo density 
profiles of DM haloes from numerical simulations is still well 
approximated by the NFW model over the range of masses 
of interest here. Motivated by this result, we therefore as- 
sume that the DM density profile is the NFW one : 



Ph(r) 



with 



AirR z s g(R v ir/R s 



(0 



i + 7 r) (is) 



g{x) = \n(l + x)-x/{l + x) . (14) 

In Eq. (|13[) . M v i r and R v i r are the virial mass and radius. 
They are not independent being related by 



Rvir 



( 3M„j r 

y-iivA th pM 



1/3 
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Figure 1. Examples of simulated rotation curves with superimposed theoretical curves. From left to right, model parameters are 
(log M d) log M vir ,c,f DM , log rix) = (11.15,12.90,10.24,0.47,0.36), (10.90,11.76, 14.77,0.45,-0.92), (10.04,12.10,13.76,0.54,1.11), while 
the simulation parameters are set as discussed in the text. Note that, depending on how the model parameters are set, it is possible to 
get rotation curves which are flat, decreasing or increasing in the outer region. 



with Ath the overdensity for spherical collapse and pu = 
3H$ n M /8irG the mean ma tter density today. We fol- 
low iBrvan fc Norman! (|l998l ) for A t h and set (Qm , h) = 
(0.28,0.70) in accordance with lKomatsu et al. I |2010T l. 

Because of the spherical symmetry, the Newtonian ro- 
tation curve may be easily evaluated as : 



GM h (r) GM vir g{r/R s ) 



Rvir Q^Rvir / Rs ) 



(15) 



While the Newtonian contributions to the rotation curve 
may be computed in the usual way, the Yukawa -like terms 
in the potential give rise to two further terms in v^(R) that 
have to be computed numerically. To this end, we must only 
insert Eqs.|TT]) and |13J) into Eq.© with / r (A ) given by 
Eq.© subtracting the first term in parentheses. Some care 
must be taken in choosing the reference radius r s . Since 
the data typically probe a limited range in Rd, a natural 
choice is to set r 3 = Rd- However, the reference density po 
is not the density at r s . It is indeed more convenient to 
set po = Pd(Rd,0) for the disc and po = Ph(Rs) for the 
halo. With such a choice, for the halo, the dimensionless 
density profile entering Eq.((9| is given by the r- dependent 
part of Ba. (|13[) provided r/R s is replaced everywhere by 
r\jr\ s . Finally, we remember the reader that, when computing 
the total rotation curve, the Newtonian terms, for both the 
disc and the halo, given by Eqs.{T5J) and (| 15p respectively, 
must be rescaled by the factor 1/(1 + 5) in order to recover 
the classical results in the GR limit (A — > oo). 



3.2 Simulating the rotation curve 

In order to be useful, our approach should rely on simulated 
rotation curves that are as realistic as possible. By this, we 
mean that a.) they must refer to spiral galaxies with reason- 
able values of the model parameters and b.) the sampling 
and the noise should be the same as actual data. Point a.) is 
the easiest to address. As a first step, we randomly generate 
the disc scalelength Rd and the halo virial mass from flat 
distributions over the ranges : 



0.5 < Ra/Rd 



< 2.0 , 11.5 < logAf„ ir s; 13.5 



with Rd.Mw = 2.55 kpc the disc scalelength for the M ilky 
Way (|Dehnen fc Binnev lll998l ; ICardone fc Sereno 1120051 ). In 
order to set the disc mass, we first define the halo mass 
fraction within the optical radius as : 



/dm 



M h (R 



opt) 



Md + Mh {R 



opt) 



(16) 



with R pt = 3.2Rd the optical radius and we have approx- 
imated the disc mass within Rd with the total disc mass. 
We then randomly generate fn m from a flat distribution in 
the range (0.9, 1. 1)/ P m. fid and /dm, fid = 50% (see, e.g., 
IWilliams et al. I (|2010t l and refs. therein). The halo scale- 
length R s is computed as R 3 = Rvir I c where the concen- 
tration c is randomly generated from a Gaussian centred 
on : 



c= 16.7 



M vi 



10 11 h- 1 M Q 



(17) 



and variance set to 10% of th e mean. Note that the a bove 
relation has been derived by iNapolitano et af~l |2005h for 
the mass r ange (0.03, 30) x 10 12 Mp) following the method 
detailed in iBullock et al. I <|200lh and updating the cosmo- 
logical model. Finally, we need to set the modified potential 
paramet ers (S, A). Following the r esult obtained for f(R)- 
gravity (|Capozziello et al~l 200911 ) . we first set 8 = 1/3 
and run different simulations randomly generating log r/\ = 
log(A/i?d) from a flat distribution covering the wide range 
(—2, 2). In order to explore the impact of 5, we also consider 
the extremal case S = 1.0 thus maximizing the contribution 
of the Yukawa terms. 

Having thus generated a realistic galaxy model, we now 
need a rule for sampling it and adding noise in such a way 
that the simulated rotation curve is similar to an observed 
one. A unique choice is not possible since the details of any 
observation depend on the instrumental setup, the observing 
conditions and the galaxy surface brightness. After a visual 
examination of rotation curves samples in literature, we have 
adopted the strategy summarized below. 

(i) Take 2Nsim equally spaced points rji in the range 
("Hmin , rj m ax) ■ and replace each of them with fji = etr/i with 
Si a randomly generated from the range (0.9, 1.1). 

(ii) For each point in the sample, generate v s im{fii) from 
a Gaussian distribution centred on the theoretical value 
and with variance set to (e c /'2)v c (fji). 

(iii) Set the error on the i-th point as ot — SiE c Vsim(jji) 
with Si randomly chosen in the range (0.9, 1.1). 
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(iv) Generate two random numbers (mi, 112) in the range 
(0, 1) and take the point i if Ui < 112- 

Typical rotation curves are well sampled up to the optical 
radius R opt = 3.2Rd, while the sampling gets worse at larger 
radii. On the contrary, the percentage errors are of the same 
order in both regions, although it is possible that the inner- 
most points have larger uncertainties because of deviations 
from ordered motions. For given input model parameters, 
we therefore generate two samples of points by first setting 

{Nsi m ,r)min,n™ax,E c ) = (30,0.1,3.1,0.25) 



and then 



CM 



= (20,3.0,10.0,0.20) 



We then add the two samples and rescale the errors in such 
a way that the standard x 2 equals 1 for the input rotation 
curve. Although such values are arbitrary, we have checked 
that the simulated rotation curves have a similar sampling 
and uncertainties of many dataset in literature (see Fig.[T|) 
so that we will not try other possible combinations. In or- 
der to quantify the bias on the model parameters, we then 
simulate 1000 rotation curves and fit different models (as 
described in the following) to determine their best fit pa- 
rameters. Note that we will not consider the errors on the 
fit quantities since they strongly depend on the uncertain- 
ties on the data points hence on the choice of the simulation 
parameters (Afsim,rimin,Vm.ax,£c)- Since we are interested 
in a statistical analysis of the full sample rather than on a 
detailed study of a particular case, we prefer to limit our 
attention to the best fit parameters only thus avoiding to 
deal with the details of the simulation procedure. 



4 BIAS ON HALO PARAMETERS 

The sample of simulated rotation curves constructed by the 
procedure detailed above is the starting point of our anal- 
ysis. Indeed, we can now fit each curve with a given (New- 
tonian) model and compare the output best fit parameters 
with the input ones. There is, however, a preliminary caveat 
to be addressed. When fitting a whatever model to a given 
dataset, one has to put down an objective criterium to deem 
the model as reliable or not. It is then common practice to 
compute the standard x 2 , i.e. 



X 



E 



Vc(Ri) - V c>t h(Ri 



and then evaluate the quality of the fit considering the value 
of x 2 = X 2 /d.o.f. with d.o.f. = Af — Af p the number of de- 
grees of freedom and Af (Afp) the number of data points (pa- 
rameters). Formally, one can then conclude that the model 
is a good description of the data if x 2 Xth with the thresh- 
old value depending on A/". Since, for our simulated curves 
Af ~ 50, Xth ~ 1-2 which is highly demanding selection 
criterium. Actually, this formal rule rigorously applies only 
if the data points are uncorrelated and the measurement 
uncertainties are Gaussian distributed. Both these assump- 
tions break down for our simulated curves since we have 
rescaled the uncertainties in such a way that X 2 = 1 for 
the input model. By a visual examination of many fitted 



curves, we have checked that a reasonably good fit (with 
no systematic deviations from the outer or inner data and 
rms percentage deviations smaller than 10%) is obtained up 
to x 2 — 2.5 which we choose as our cut. We will therefore 
consider a model as successfully fitting the data if x 2 ^ 2.5 
and log M V i r ^ 14 with this latter criterium introduced to 
avoid unphysical solutions. However, we will also discuss the 
trends for the quantities of interest for a subsample obtained 
by using a stringent criterium, i.e. x 2 ^ 1.5 in order to ex- 
plore the impact of the threshold x 2 value. 



4.1 Navarro-Frenk- White models 

In realistic situations, one has a set of (7?, v c ) data and a 
measurement of the galaxy surface brightness in a given 
band. It is then common to set the disc scalelength Rd to the 
value inferred from photometry, while the disc mass can be 
inferred from the total luminosity provided an estimate of 
the stellar M/L ratio is somewhat available (e.g., from stel- 
lar population synthesis models fitted to the galaxy colours). 
As a first step, we therefore assume that both (Rd, Md) are 
known and fit the simulated data for each rotation curve to 
determine the NFW model parameter^ only. 



4-1.1 5 = 1/3 simulated samples 

Let us first consider the results obtained setting 5=1/3 
when simulating the rotation curves. As yet said above, we 
are here trying to fit modified gravity rotation curves us- 
ing theoretical models computed in a Newtonian framework 
so that the first point to address is whether this is indeed 
possible. To this end, it is sufficient to check what is the per- 
centage of rotation curves passing the two selection criteria 
quoted before. We indeed find that 90% of the simulated 
rotation curves are fitted with x 2 ^ 2.5 and logAf V f r ^ 14 
(which we will refer to as the well fitted sample, hereafter 
WS), while this fraction decreases to 77% if x 2 ^ 1.5 is de- 
manded (defining what we will hereafter refer to as the best 
fitted sample, BS). Averaging over the sample values gives : 

(x) = 1-27 ± 0.24 (1.19 ± 0.13) for WS (BS) 

while for the rms of the percentage deviations Av c /v c = 
(v° ba - vl u )/v° c ba we get : 

(rms(Av c /v c )) = 6.4% ± 0.8% (6.3% ± 0.7%) for WS (BS) . 

Motivated by the high fraction of successful fits and the 
low values of x 2 an d rms(Av c /v c ), we can therefore safely 
conclude that it is possible to fit a modified gravity rota- 
tion curve with a NFW Newtonian one obtaining both an 
acceptable x 2 value and reasonable model parameters. 

Having checked that the model reasonably well fit the 
data, we can rely on the best fit model parameters and com- 
pare them with the input ones to estimate what is the bias 



5 To this end, we use a straightforward Monte Carlo Markov 
Chain algorithm to minimize the \ 2 with respect to the param- 
eters (log 77s, log V v i r ) and then infer the values of (c,M v i r ) and 
the dark matter mass fraction foM- This approach is more stable 
than fitting directly for (c, M v i r ) and allows to correctly recover 
the input model parameters when we fit rotation curves simulated 
in the Newtonian framework. 
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Figure 2. lZ(x) vs log 77^ with x = r] s (upper left), M v i r (upper right), c (lower left), foM (lower right) for the NFW model fit to the 
simulated rotation curves with 5 = 1/3. Note that only 10% of the WS sample is plotted to not clutter the figure. 



Table 1. Bias 1Z(x) on the NFW model parameters assuming the disc mass is known and 5=1/3. Columns are as follows : 1. parameter 
id; 2., 3., 4. mean ± standard deviation, median and rms values, 5., 6., 7., 8., 9. Spearman rank correlation coefficients between TZ(x) 
and input logjjs, logM„i r , c, fnMt log^A- Upper half of the table is for the WS sample, while lower half for the BS one. 



Id 


(K) 


Turned 




c(iog Vs ,n) 


C(\ogM vir ,K) 


C{c,K) 


c(f DM ,n) 


c(io gVx ,n) 




4.0 ±2.8 


2.9 


4.8 


0.31 


0.18 


-0.17 


-0.09 


-0.47 




3.9 ±3.8 


2.4 


5.5 


0.48 


0.31 


-0.27 


-0.05 


-0.32 


C 


0.45 ±0.15 


0.45 


0.48 


-0.13 


-0.05 


0.07 


0.13 


0.55 


Idm 


0.70 ±0.06 


0.70 


0.71 


0.36 


0.30 


-0.20 


0.39 


0.70 


Is 


3.8 ±2.8 


2.8 


4.7 


0.32 


0.24 


-0.23 


-0.07 


-0.42 




3.6 ±3.6 


2.2 


5.1 


0.50 


0.36 


-0.34 


-0.02 


-0.26 


C 


0.47 ±0.15 


0.47 


0.49 


-0.14 


-0.10 


0.13 


0.12 


0.53 


Jdm 


0.70 ±0.06 


0.70 


0.71 


0.39 


0.31 


-0.18 


0.40 


0.71 



induced by an incorrect assumption of the gravity law. To 
quantify this bias, we define IZ(x) = x fu/ 'x s im, i.e. the ratio 
between the best fit and the input values of a given quan- 
tity x. Needless to say, should IZ(x) be on average close to 
1, we can conclude that assuming Newtonian gravity to fit 
modified gravity rotation curves does not induce a signifi- 
cant bias on the estimate of x. The results, summarized in 
Table [T] for both the WS and BS samples, show that such 
a bias is indeed quite important and only mildly depend on 
the selection criteria adopted (so that we will hereafter refer 
to the WS results only unless otherwise stated). 

Table[T]gives some statistics on the distribution of lZ(x) 
for the halo parameters (rj s , M v i r ). However, these values 



could be misleading since the histograms for IZljis) and 

lZ(M V i r ) are strongly asymmetric with long tails towards 

the right. In other words, TZ{r) s ) and lZ(M v i r ) could be much 

larger than their median value with lZ(r] s ) being as high 

as 12 and values of lZ(M V i r ) as large as 15 so that both 

the halo scalelength and virial mass can be grossly over- 

1/3 

estimated. Since c oc R V i r oc M vir , one could naively ex- 
pect that the concentration is overestimated too. On the 
contrary, the distribution of the lZ(c) values is almost sym- 
metric around its mean clearly disfavouring IZ(c) > 1, i.e. 
the concentration is underestimated. Actually, such a re- 
sult can be understood remembering that c oc R^ 1 so that 
11(c) oc 1Z 1 / 3 (M vir ) /1Z(rj s ) finally leading to values smaller 
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than unity. The emerging picture is therefore that of a halo 
having a larger mass than the input one, but also a larger 
scalelength. This can be qualitatively explained noting that 
the Yukawa -terms in the rotation curve increases the net 
circular velocity both in the inner and outer regions probed 
by the data. In order to adjust the fit in the halo dominated 
regions, one has to increase M v i r , but then R s has to be 
increased to in order to lower the density (and hence the 
contribution to the rotation curve) in the inner regions not 
to overcome the observed circular velocity. As a result, the 
concentration is underestimated and the same takes place 
for the dark matter mass fraction within the optical radius 
since the halo mass is now pushed outside the optical ra- 
dius thus explaining why 1Z(}dm) is smaller than unity. As 
Fig-El shows, IZ(x) is correlated with log??}, with the sign of 
the correlation depending on lZ{x) being typically larger or 
smaller than 1. Not surprisingly, the larger is the Yukawa 
scalelength A, the closer is IZ(x) to 1, i.e., the smaller is the 
bias induced by the assumption of Newtonian gravity. How- 
ever, there is a large scatter in these correlations since the 
relative importance of the correction also depends on the 
input halo parameters. Indeed, when R 3 ~ A, the halo con- 
tribute to the modified rotation curve is maximized so that 
the Yukawa corrections are more important and the bias of 
the fitted parameters is larger. As a consequence, one can 
therefore expect a correlation between IZ(x) and log rj s which 
is what we indeed find looking at the results in Table \T\ 

4-1-2 S = 1.0 simulated samples 

Let us now consider the results for the rotation curves sam- 
ple simulated under the extreme assumption 5 = 1.0, i.e. 
when the contribution of the Yukawa term to the point mass 
potential Q is the same order of magnitude as the Newto- 
nian one. Note that, different from the previous case, S = 1.0 
is, as far as we know it, an unmotivated assumption, but we 
consider it to investigate the impact of 5 on the results. 

The first striking result is that the fraction of well fit- 
ted rotation curves dramatically drops to 22% applying the 
WS selection criteria (x 2 ^ 2.5 and logM„i r ^ 14), while no 
curves pass the cut x 2 ^ 1-5 thus showing that it is likely 
not possible to reproduce the simulated curves with NFW 
model if Newtonian gravity is assumed. Such a result can 
be understood noting that, for 5 = 1.0, the modified grav- 
ity rotation curve is much larger than the Newtonian one 
so that one has to greatly increase the halo virial mass to 
fit the outer rotation curve. In order to compensate the cor- 
responding increase in the inner region, one should set r/ s 
as large as possible, but this also tend to lower the circu- 
lar velocity for R » R a . Finding a compromise between 
these two opposite trends is quite difficult thus leading to 
unacceptably large x 2 values. 

It is worth investigating what is the bias induced on 
the halo model parameters for the successfully fitted cases. 
The trends of the bias parameters lZ(x) are the same as 
in the previous case, but the typical values are now much 
larger. In particular, TZ(-q s ) can be as high as 40 thus leading 
to grossly underestimates of both the concentration c and 
the dark matter mass fraction /dm- Not surprisingly, the 
only way to reduce such large biases is to increase log^A 
to values larger than our upper limit log r]\ = 2. However, 
in this case, the modified gravity potential differs from the 



Newtonian one only on group and cluster scales where the 
rotation curve are no more measured so that our analysis 
becomes meaningless. 



4.2 Generalized Navarro-Frenk- White models 

Up to now, we have considered the NFW profile for the 
halo model, but such a choice does not allow us to inves- 
tigate whether the mismatch between modified and Newto- 
nian gravity can have an impact on the inner logarithmic 
slope. In order to explore this issue, we therefore consider 
the generalized NFW model (gNFW, Jing & Suto 2000) : 



Ph(r) = 



with 



M vi 



h(x) 



-InRlhiRvir/Rs) 



1 + i 



-(3-c) 



(i + C) 3 



(18) 



(19) 



The parameter a depends on the behaviour of the rotation 
curve in the inner regions where the disc contribute can be 
larger than the halo one. It is therefore worth investigating 
which is the impact of the uncertainties on Md in order to 
see how the bias on a depend on it. 



4-2.1 gNFW models with known disc mass 

As a first step, we assume that Md is set to the input value 
and only look at the bias on the parameters (a, c, M„; r ). As 
a preliminary remark, it is mandatory explaining how the 
concentration is defined for gNFW models. To this end, one 
may simply note that, for the NFW density profile, R s is the 
radius at which the logarithmic slope equals the isothermal 
value 7 = —2. Therefore, the concentration for the gNFW 
model may be easily defined as : 

Rvir Rvir (20) 



R- 



(2 - a)R a 



with 7(-R_ 2 ) = -2. Note that, for a = 1, the gNFW model 
reduces to the NFW one and Eq. (|20|) gives the standard 
definition for the concentration of NFW haloes. 

Applying the selection criteria used above to the sample 
with S = 1/3, we find that the WS sample contains now 90% 
of the simulated curves, while this fraction lowers to 82% for 
the BS sample. Averaging over the samples gives : 

(X 2 ) = 1.28 ± 0.19 (1.23 ± 0.12) for WS (BS) , 

(rms(Av c /v c )) = 6.4% ± 0.8% (6.3% ± 0.8%) for WS (BS) . 

These values are almost identical to those obtained for the 
fits of the NFW model described before and allows us to 
confidently conclude that it is possible to reproduce the sim- 
ulated rotation curves with a gNFW model and Newtonian 
gravity. In a sense, this is not surprising since the gNFW 
model has one more parameter than the NFW one so that 
there is a much larger freedom to adjust the shape of the 
theoretical rotation curve to fit the observed one. As a re- 
sult, the percentage of well fitted rotation curves slightly 
increases, but the quality estimator (x 2 and the rms value 
of Av c /v c ) stay almost unchanged. 

The bias induced on the model parameters can be still 
estimated considering the quantity IZ(x) defined before and 
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Figure 3. Histogram of 1Z(a) values and its dependence on logr;>, for the gNFW model (with disc mass set to the input value) fit to 
the simulated rotation curves with <5 = 1/3. Note that only 10% of the WS sample is plotted to not clutter the figure. 



Table 2. As in Table [T] for the gNFW with fixed disc mass case. 



id 


(n) 


Turned 




C{logVs,n) 


C(\ogM vir ,Tl) 


C(c, 71) 


C(fDM,K) 


C(log Vx ,n) 


a 


0.92 ±0.22 


0.95 


0.95 


-0.56 


-0.58 


0.37 


0.09 


0.32 


Is 


2.8 ±2.6 


2.0 


3.8 


-0.11 


-0.06 


0.0 


-0.05 


-0.06 




2.4 ± 2.6 


1.4 


3.5 


0.43 


0.39 


-0.39 


0.03 


-0.18 


C 


0.56 ±0.18 


0.53 


0.59 


0.15 


0.08 


-0.03 


-0.03 


0.16 


Idm 


0.70 ±0.06 


0.71 


0.71 


-0.06 


-0.18 


0.04 


0.40 


0.70 


a 


0.93 ±0.22 


0.96 


0.95 


-0.56 


-0.50 


0.39 


0.14 


0.31 


Vs 


2.8 ±2.5 


2.1 


3.7 


-0.08 


-0.07 


0.0 


0.11 


-0.07 




2.4 ± 2.6 


1.4 


2.7 


0.42 


0.36 


-0.37 


0.05 


-0.17 


C 


0.55 ±0.18 


0.53 


0.58 


0.12 


0.08 


-0.03 


-0.09 


0.18 


Jdm 


0.70 ±0.06 


0.70 


0.70 


-0.05 


0.18 


0.03 


0.42 


0.69 



summarized in Table [2] for the present castjf]- Comparing to 
the results in Table [T] it is apparent that the bias induced 
on the halo parameters in common, i.e. (rj 3 ,c, M V i r , /dm), 
is almost the same both for the WS and BS samples. The 
same qualitative discussion can be repeated here so that we 
will not consider anymore this issue. It is, on the contrary, 
more interesting to look at the bias on a which asks for some 
caution. As can be seen from the left panel in Fig. [3] the dis- 
tribution of a values is quite large and asymmetric, while the 
right panel is a clear evidence that its mean value strongly 
depend on log?7,\. Indeed, should we have only fitted models 
with A ^ Rd, Ti{ct) would have been smaller than 1, i.e. the 
inner slope would have been underestimated and shallower 
models be preferred. Should the modified gravity scalelength 
be smaller than the typical disc one, fitting rotation curves 
assuming the validity of Newtonian gravity would have us 
led incorrectly to believe that the inner slope of the density 
profile is much smaller than the NFW one. Such a result, 
therefore, suggests that the cusp/core controversy is just a 
fake problem originated by an incorrect assumption of what 
the underlying gravity theory actually is. 

As a final remark, it is worth noting that a second dif- 
ference among the results in Tables [T] and [2] is represented 



6 Note that, for all the simulated curves, it is a = 1 since we 
have assumed a NFW model in the simulation process. Therefore, 
11(a) = a for the gNFW models. 



by the markedly different values of the Spearman correla- 
tion coefficients. Actually, this is a consequence of having 
added one more parameter which introduces a degeneracy 
between the fitted quantities hence partially washing out the 
correlations of IZ(x) with the input halo parameters and the 
modified gravity scalelength. In particular, a and logr/ s are 
clearly correlated since both set the shape of the inner rota- 
tion curve. Although with a large scatter, IZ(a) is typically 
smaller (larger) than 1 for log?7A < (> 0) so that 1Z(at) 
actually correlates with log^A, but not linearly thus giving 
rise to a small Spearman correlation coefficient (which looks 
for linear correlations). Since r] a has to be adjusted accord- 
ing to the changes in a, lZ(r] a ) has to follow lZ(a) thus losing 
its correlation with log77,\. 

Such results are strengthened if we consider the simu- 
lated curves for S = 1.0, but now relying on a much smaller 
statistics. Indeed, the WS sample is now made out of only 
20% of the full set of simulated curves, while only 2% of 
them are included in the BS sample. These fractions are 
larger than in the NFW fits considered before, but are still so 
small to make us conclude that such extreme modified grav- 
ity curves can not be reproduced in a Newtonian framework 
using the gNFW model. For the few surviving curves, we 
can repeat the same discussion above with the only caveat 
that now the IZ(x) values deviate more from 1 than in the 
5=1/3 case. Moreover, the slope a can be underestimated 
also when log?7A > depending on the input halo parame- 
ters. This is a still further evidence in favour of the suggested 
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Figure 4. Histograms of 1Z(x) with x = (left) and x = a (right) for the gNFW model with free disc mass fits and <5 = 1/3. 



Table 3. As in Table [2] for the gNFW with free disc mass case. 



Id 


(11) 






C (log M d ,1Z) 


C(logr, s ,7^) 


c(io g M vir ,n) 


C(c, K) 


C(f DM ,K) 


C (log Vx , 11) 


M d 


0.79 ±0.19 


0.80 


0.82 


0.0 


0.14 


0.12 


-0.07 


-0.02 


0.09 


a 


1.00 ±0.20 


1.02 


1.01 


-0.10 


-0.01 


-0.08 


-0.01 


0.01 


0.33 


Is 


1.6 ± 1.5 


1.2 


1.6 


-0.01 


0.21 


0.18 


-0.15 


0.03 


0.0 




1.7 ±3.5 


0.8 


3.8 


0.06 


0.22 


0.23 


-0.18 


0.02 


-0.09 


C 


0.87 ±0.37 


0.82 


0.95 


0.0 


-0.22 


-0.21 


0.16 


-0.03 


0.07 


}dm 


0.97 ±0.24 


0.96 


1.00 


-0.04 


-0.12 


-0.12 


0.07 


0.05 


0.10 


M d 


0.78 ±0.20 


0.78 


0.80 


0.03 


0.11 


0.13 


-0.09 


-0.04 


0.08 


a 


1.00 ±0.20 


1.02 


1.02 


-0.12 


0.03 


-0.06 


-0.02 


0.02 


0.35 


Is 


1.6 ± 1.6 


1.2 


2.2 


0.00 


0.20 


0.18 


-0.15 


0.01 


-0.01 




1.6 ±3.2 


0.8 


3.5 


0.07 


0.21 


0.23 


-0.18 


0.0 


-0.12 


C 


0.89 ±0.38 


0.84 


0.97 


-0.02 


-0.22 


-0.22 


0.17 


-0.01 


0.09 


Jdm 


0.99 ± 0.24 


0.97 


1.01 


-0.06 


-0.09 


-0.12 


0.08 


0.07 


0.10 



interpretation of the cusp/core controversy as the outcome 
of a systematic error in the assumed gravity theory. 

4-2.2 gNFW models with free disc mass 

Up to now, we have set the disc mass to the input value im- 
plicitly assuming that one is able to perfectly estimate this 
quantity from the galaxy colors and a stellar population syn- 
thesis code. Actually, this is not always the case and, on the 
contrary, it is usual to include M d in the set of parameters 
to be determined by fitting the rotation curve data. We have 
therefore repeated the above analysis for the gNFW model 
leaving M d free to be adjusted by the fit. 

Considering there is one more fit parameter, it is not 
surprising that the fraction of successfully fitted rotation 
curves (for 5=1/3) increases to 92% (85&) for the WS (BS) 
samples and also the quality of the fit is improved being : 

(x) = 1.19 ± 0.20 (1.14 ± 0.13) for WS (BS) , 

(rms(Av c /v c )) = 6.1% ± 0.6% (6.0% ± 0.8%) for WS (BS) . 

Table [3] summarizes the results on the bias estimator TZ(x) 
for the different parameters involved which allow us to draw 
some interesting considerations. First, we note that the halo 
model parameters are now better recovered than in previous 
cases. For instance, although the corresponding distributions 
have long tails up to IZ(x) ~ 10 — 15, most of the values of 



1Z(r] s ) and TZ(M v i r ) are smaller than 2 leading to a modal 
value (i.e., the most peak of the histogram) close to 1, i.e. 
{rj s ,M v i r ) axe not biased anymore for most of the WS and 
BS sample curves. A similar discussion also apply to the con- 
centration c and the dark matter mass fraction /dm whose 
histograms are quite symmetric and centred close to the no 
bias value, 1Z(x) = 1. Concerning a, one must still keep in 
mind that its mean value depend on the cut on log rj\ typi- 
cally being 1Z(a) < 1 (i.e., the inner slope is shallower than 
the input one) when A < R d . Note, however, that this effect 
is now less pronounced than before thus explaining why the 
distribution of TZ(a) values is peaked close to 1 with a not 
too large width. On the contrary, the disc mass M d turns 
out to be biased low with M d being smaller than the input 
value for most of the cases. Such an unexpected result can be 
qualitatively understood considering that the Yukawa terms 
make the modified gravity rotation curve larger than the 
Newtonian one in the inner and outer regions. As in the 
case of the NFW fit, one must increase the halo virial mass 
to recover this additional contribution in the halo dominated 
regions which increases also the halo inner circular velocity. 
In the NFW case, one has then to increase r) s to prevent 
overestimating the simulated curve, while the same effect is 
now accomplished by lowering the disc mass. This explains 
why the halo parameters turn out to be less biased than in 
the NFW fit case. It is, however, worth stressing that, as 
shown in the left panel of Fig. [4] the lZ(M d ) distribution is 
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actually bimodal. As a consequence, when the disc mass is 
not biased low, one has again to increase 77,, thus explaining 
the long tails towards TZ(rj s ) » 1 values. Alternatively, one 
can leave rjs unchanged, but make the density profile shal- 
lower hence giving rise to the tail towards small lZ(a) < 1 
in the right panel of Fig.0] Model degeneracies also explain 
why the correlation coefficients among IZ(x) and the input 
model parameters are typically quite small. We, however, 
caution the reader that a small value of Spearman correla- 
tion does not necessarily imply that IZ(x) doe not depend 
on the corresponding quantity (see, e.g., the IZ(a) depen- 
dence on log77A for the gNFW fits with fixed disc mass). 
We have therefore checked whether this is the case or not 
by directly looking at the 1Z(x) vs pi plots (with pi the i- 
th input parameter) finding that the scatter of the points in 
each plot clearly indicates a complicated dependence a weak 
dependence on all the input parameters but log 77a. 

Finally, we have repeated the above analysis setting 
5 = 1.0 to generate the simulated rotation curves. Differ- 
ently from the previous cases considered up to now, we now 
find that it is possible to fit the gNFW model with free disc 
mass to the simulated rotation curves. In a sense, the addi- 
tional freedom we have now to adjust the disc mass makes 
it possible to recover the inner rotation curve better than in 
the previous models thus leading to a successful fit for 97% 
(85%) of the galaxies using the selection criteria for sample 
WS (BS). The quality of the fit is also remarkably good : 

(X 2 ) = 1.32 ± 0.30 (1.18 ± 0.15) for WS (BS) , 

(rms(Av c /v c )) = 6.5% ± 1.0% (6.2% ± 0.7%) for WS (BS) . 
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As it is clear from the large standard deviations and rms 
values, the distribution oflZ(Md) and lZ(M V i r ) are strongly 
asymmetric and, as a result, the one for a presents a non 
negligible amount of points with TZ(a) > 1. This result can 
be explained noting that setting 5 = 1.0 maximizes the con- 
tribution of the Yukawa terms so that we now have to mimic 
a similar effect by adjusting the model parameters in the 
Newtonian fitting. One possibility is to leave a unchanged 
with respect to the input value (a = 1.0), but increasing 
the halo mass by a significant factor thus leading to large 
values of lZ(M V i r ). But, in this case, one has also to decrease 
the disc contribution to the inner rotation curve to compen- 
sate for the additional dark matter. On the contrary, one 
can leave almost unchanged M V i r , but redistribute the dark 
mass pushing it towards the inner region. To this aim, one 
should make the profile steeper, i.e. increase a, while Md 



must be smaller to not overcome the inner circular veloc- 
ity. This strategy will lead to IZ(a) > 1 and 1Z(M V i r ~ 1, 
while again it is IZ(Md) < 1. The final histograms of lZ(x) 
for {Md, a, M v i r ) is then the outcome of a mixture of these 
two possible strategies whose relative contribution depends 
on the details of the individual rotation curves. 



5 BIAS ON THE HALO DENSITY PROFILE 

The use of the gNFW model to fit the simulated rotation 
curve may also be read as a first step towards completely 
mismatching the halo density profile. In a sense, the gNFW 
model departs from the input NFW one only because of 
the possibly different inner logarithmic slope. It is, however, 
common in data analysis to use also models that differs from 
the NFW one both in the inner and outer regions. Moreover, 
the gNFW is still a cuspy model, while the cusp/core contro- 
versy comes from the observation that cored models better 
fit the observed rotation curves. To this end, we now drop 
the implicit assumption that the trial model tracks or gener- 
alizes the NFW profile and investigate whether completely 
different models in the Newtonian framework are in accor- 
dance with our modified gravity simulated rotation curves. 



5.1 The pseudo - isothermal sphere 

A classical example of a cored model is represented by the 
pseudo - isothermal (hereafter, PI) model whose density pro- 
file reads : 



Concerning the values of TZ(x) quantities, they are almost 

the same as those in Table[3]for the parameters (773 , c, /dm) Ph( r ) = 

which are therefore almost unbiased. On the contrary, for 

the other parameters (and the WS sample), we get : with 



M vi 



^■KR 3 a h pi (R vir /R s 



1 + 



r 2 
R~s 



h(x) = x — arctan x 



(21) 



(22) 



As it is clear, the PI density law approaches a constant value 
for r « R s , while falls off as r~ 2 in the outer regions. As 
such, the PI model is radically different from the NFW and 
gNFW ones so that it is interesting to see whether can fit or 
not the data. Following common practice, we leave the disc 
mass as an unknown and estimate (Md,R s , M v i r ). 

Applying the corresponding selection cuts, we find that 
78% of the simulated galaxies fall into the WS sample, while 
the BS one contains 43% of the full sample. Such high frac- 
tion and the low values of both the reduced \ 2 

(X 2 ) = 1.49 ± 0.35 (1.24 ± 0.15) for WS (BS) , 

and of the rms of percentage residuals 

(rms(Av c /v c )) = 7.0% ± 0.9% (6.9% ± 0.9%) for WS (BS) 

make us confident that the Newtonian circular velocity of 
the PI model can indeed mimic the rotation curve predicted 
by modified gravity. This is in agreement with what is indeed 
found in the literature where galaxies rotation curves data 
are typically well fitted by PI haloes. 

Having used now a different halo model to fit the data, 
it is not possible to straightforwardly compare the output 
parameters with the input ones. First, although we use the 
same symbol, the scalelength R s of the PI model is not de- 
fined by the condition 'y(Rs) = —2 as for the NFW model, 



12 V.F. Cardone & S. Capozziello 



Table 4. As in Table [T] for the PI model with free disc mass case. 



Id 


(n) 


Turned 




C (log M d ,K) 


C(logr, s ,ft) 


C '(log M vir ,K) 


C(c, K) 


C(f DM ,K) 


C(log7? A ,7£) 


M d 


1.20 ±0.10 


1.20 


1.20 


-0.12 


-0.26 


-0.21 


0.09 


0.55 


0.48 




5.2 ±7.0 


3.8 


8.7 


0.02 


0.12 


0.09 


-0.03 


0.0 


0.06 




0.37 ±0.08 


0.38 


0.38 


0.08 


0.49 


0.38 


-0.30 


-0.06 


-0.22 


M d 


1.19 ±0.10 


1.19 


1.19 


-0.13 


-0.27 


-0.19 


0.15 


0.57 


0.44 




5.4 ±8.9 


3.7 


10.0 


0.0 


0.06 


0.04 


0.0 


-0.03 


0.09 


Sum 


0.38 ± 0.08 


0.39 


0.39 


0.07 


0.56 


0.42 


-0.34 


-0.10 


-0.27 



Table 5. As inTablc [T] for the Burkort model with free disc mass case. 
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(K) 


"Turned 




C(\ogM d ,K) 


C(log Vs ,K) 


C '(log M vir , 11) 


C(c,K) 


C(f DM ,K) 


C(\o gV x,n) 


M d 


1.14 ±0.11 


1.13 


1.15 


0.0 


-0.09 


-0.03 


-0.02 


0.32 


0.30 




1.1 ± 1.1 


0.7 


0.8 


0.06 


0.13 


0.16 


-0.15 


0.09 


0.13 


}dm 


0.46 ±0.11 


0.44 


0.48 


-0.06 


-0.03 


-0.07 


0.11 


0.05 


-0.08 


M d 


1.11 ±0.10 


1.10 


1.12 


-0.11 


0.03 


0.05 


0.05 


0.33 


0.49 




1.0 ±1.1 


0.6 


1.5 


0.05 


0.12 


0.14 


-0.10 


-0.10 


0.21 




0.48 ±0.11 


0.47 


0.50 


-0.01 


-0.05 


-0.07 


-0.01 


0.11 


-0.22 



but rather represents the core radius. As such, it is meaning- 
less to compare the input n s with the output R s /R d since 
they refer to a different characteristic of the density profile. 
Moreover, a concentration for the PI model is not defined 
so that we can not compare with the input one. On the con- 
trary, the total disc mass M d , the virial mass M v i r and the 
dark matter mass fraction fr>M refer to global properties 
of the disc + halo model so that it makes sense to compare 
them with the input ones. Table|4]then gives the bias on the 
global parameters (M d , M v i r , /dm) both the WS and BS 
samples (having set 8 = 1/3). It is clear that all these three 
quantities are severely biased with (M d , M v i r ) being overes- 
timated and /dm grossly underestimated. Such a behaviour 
can be explained noting that, once the core radius has been 
set, the only way to increase the Newtonian rotation curve 
is to increase the global masses (M d , M V i r ) thus leading to 
asymmetric distributions both peaked in 1Z(x) > 1. It is 
worth noting that, if one ignores that the underlying grav- 
ity theory is not Newtonian, a large disc mass will be in- 
terpreted as the presence of a maxi mal disc which is in- 
deed the case in many spiral galaxies jPalunas fc Williams I 
l200d ). Being both M V i r biased high, one could find the re- 
sult IZ(fDAi) << 1 an unexpected nonsense. Actually, one 
has first to note that also the disc mass is overestimated 
which automatically reduces the dark matter mass fraction. 
Moreover, for the PI model, the mass profile increases al- 
most linearly, while, for the input NFW profile, Mh(r) grows 
approximately in a logarithmic way. Being the input R 3 al- 
most the same as the output core radius, we therefore get 
MNFw(Ropt) > Mpi(Ro P t) even if the virial mass of the PI 
model is larger than the NFW input one. As a consequence, 
fr>M turns out to be underestimated as we indeed find. 

As a final test, let us consider what happens when 8 = 
1.0 is used in simulating the rotation curves. It turns out that 
the percentage of curves entering the WS sample reduces to 
67%, while only a modest 12% enter the BS sample. This 
sudden drops is actually due to the cut on the output virial 



mass rather than on the reduced \ 2 ■ As such, we believe 
that such a case can not be mimicked by a PI model under 
the Newtonian gravity assumption and not discuss anymore 
the bias on the output parameters. 

5.2 The Burkert model 

Another cored model but with a different o uter profile is t he 
Burkert model whose density profile read (|Burkert Ill995t ) : 

Ph(r)= 4,mh^R mr /R s ) { 1+ Rz) \ 1 + i) 2 (23) 
with 

h B {x) = ln(l + a:) - arctana;+ (1/2) ln(l + x 2 ) . (24) 

Note that the Burkert model presents an inner core with ra- 
dius R s , but asymptotically drops off as r -3 (hence having 
a finite total mass). As such, it represents a sort of compro- 
mise between the cored PI model and the cusped NFW one. 
Again, we will leave the disc mass free and determine the 
three parameters (M d , r] a , M v i r ) from the fit. 

Setting 8 = 1/3, we find that 88% of the galaxies fall 
into the WS sample, while this fraction drops to 37% for 
the BS one. The quality of the fit may be judged from the 
following figures of merit : 

(x) = 1-63 ± 0.32 (1.31 ± 0.15) for WS (BS) , 

(rms(Av c /v c )) = 7.3% ± 1.0% (6.4% ± 0.7%) for WS (BS) , 

comparable to the values obtained for the PI fits. It is worth 
noting that the fraction of galaxies in the WS sample is 
larger than for the PI model as a result of the different mass 
profile. Indeed, for the PI model, asymptotically M(r) oc r 
so that v c (r) oc M(r)/r ~ const, while for the Burkert model 
we have a finite total mass and hence a Keplerian fall off. As 
Fig. [1] shows, our simulated rotation curves sample is made 
out of galaxies which can be flat, decreasing or increasing in 
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the outer regions depending on the input model parameters. 
It is, of course, quite difficult to fit decreasing curves with the 
PI model which does not predict such a behaviour, although 
the limited radial range probed and the uncertainties allow 
to get a reasonable match in some cases. Models with ph oc 
r -3 for r >> R s work better thus motivating the higher 
fraction of success of the Burkert profile. 

The bias values lZ(x) axe summarized in Table[5]for the 
case with 5 = 1/3 and shows that, even assuming a Burk- 
ert halo, the disc mass turns out to be overestimated hence 
mimicking maximal disc solutions. On the contrary, the halo 
virial mass is only modestly biased, especially if compared 
to the outcome of previous fits. Although the lZ(M v i r ) dis- 
tribution has a long tail to the right of the mean value, 
TZ(AI V ir) ^ 2 for most of the galaxies in the WS sample. 
Since /d m depends on 1 /Md and Md is overestimated (while 
the dependence on Mh(R op t) is weaker being this quantity 
present both at the numerator and denominator), it is then 
expected that Jdm is biased low. This is indeed what we 
find in accordance with the similar result obtained for the 
PI model. Note, however, that this time the bias is smaller 
than for the PI model, even if still quite significant. 

Finally, we have repeated the above analysis setting S = 
1.0 when simulating the rotation curves. The WS sample 
still contains 73% of the full simulated sample, but only 
17% of them pass the criteria for entering the BS sample. 
This situation is quite similar to what takes place for the PI 
model so that we do not discuss it here anymore. 



6 CONCLUSIONS 

General Relativity has been experimentally tested on scales 
up to the Solar System one so that assuming it still holds 
on much larger scales, such as the galactic and cosmological 
ones, is actually nothing else but an extrapolation. Moti- 
vated by this consideration and the difficulties in explaining 
the observed accelerated expansion without introducing new 
unknown ingredients like dark energy, a great interest has 
been recently devoted to modified gravity theories which 
have proven to work remarkably well in fitting the data and 
predicting the correct growth of structures. Modifying Gen- 
eral Relativity has impact at all scales so that, provided no 
departures from standard results, well established at Solar 
System scales, one cannot exclude a priori that the gravita- 
tional potential generated by a point mass source has not 
the usual Keplerian fall off, oc 1/r, but a weaker one. Here 
we have considered the case of a Yukawa- like correction, 
i.e. (j> oc (l/r)[l + 5exp(— r/A)] where the scale length A is 
related to the effective scalar field (coupled with matter) in- 
troduced by several modified theory of gravity. In particular, 
this kind of potential comes out in the weak field limit of 
/(ii)-gravity. Provided A is much larger than the Solar Sys- 
tem scale, the corrections to the potential can significantly 
boost the circular velocity for an extended system like a 
spiral galaxy. Assuming Newtonian gravity to compute the 
theoretical rotation curve to the observed one then intro- 
duces a systematic error which can bias the estimate of the 
galaxy parameters thus leading to misleading conclusions. 

To investigate this issue, we have fitted the Newtonian 
circular velocity of some widely used halo models to a large 
sample of simulated rotation curves computed using the as- 



sumed modified potential and a disc + NFW profile. Com- 
paring the input parameters with the best fit ones allows us 
to draw some interesting lessons on the consequences of a 
systematic error in the adopted gravity theory. As a general 
result, we find that, for cusped halo models, the disc mass is 
underestimated, while the halo scalelength and virial mass 
are biased high. Since the concentration c is also biased low, 
the c- M v i r relation will have a different slope and intercept 
than the one we have used to generate the model based on 
the outcome of N- body simulations. Moreover, if left free in 
the fit, the inner slope of the density profile turns out to be 
biased low if A is smaller than the disc half mass radius Rd- 
It is interesting to note that halo models shallower than the 
NFW one in the inner regions with (c, M v i r ) values not con- 
sistent with the prediction of N - body simulations and with 
maximal discs are a common outcome of fitting observed 
(not simulated) rotation curves. Such results are usually in- 
terpreted as failures of the ACDM model due to misleading 
assumptions on the dark matter particles properties (such 
as their being hot or cold and their interaction cross sec- 
tions) or to having neglected the physics of baryons in the 
simulations. Our analysis point towards a different expla- 
nation considering these inconsistencies as the outcome of 
forcing the gravitational potential to be Newtonian when 
it is not. In order to test such an hypothesis, one should 
fit a homogenous set of well sampled and radially extended 
rotation curves assuming a theoretically motivated modi- 
fied gravity potential (to set the strength 8 of the corrective 
term) and a classical NFW model (since it is in agreement 
with N-body simulations also in fourth order theories). 

Should the scalelength A of the modified potential be 
smaller than the disc one -Rd, we have shown that forcing 
the potential to be Newtonian may also lead to completely 
mismatch not only the model parameters, but also the halo 
density profiles. Indeed, cored models, such as the PI and 
Burkert ones, turn out to be able to fit equally well the 
rotation curve data. We have not carried out here a case- 
by-case comparison among the different models considered 
since this will depend critically on what the actual (not the 
simulated) uncertainties are and on the sampling and ra- 
dial extent of the data. We can however anticipate that, 
for most cases, the cored models will work better than the 
cusped ones since they are better able to increase the in- 
ner rotation curve redistributing the total dark matter mass 
inside the core thus leaving almost unchanged the outer ro- 
tation curve. Such a result may have deep implications on 
the cusp/core controversy which should then be read as an 
evidence of an inconsistent assumption about the correct un- 
derlying gravity theory. From an observational point of view, 
one could try to fit our NFW + modified potential to LSB 
galaxies since this dark matter dominated systems are usu- 
ally considered the best examples of the cusp/core problem. 
Note that our simulated sample does not actually contain 
LSB - like galaxies since we have set the input dark matter 
mass fraction as /dm ~ 50% and adjusted the disc parame- 
ters having in mind a Milky Way - like spiral galaxy. Should 
we have simulated only LSB -like systems, we expect that 
cored models will definitely be preferred over cusped ones 
thus further strengthening our conclusions. 

As a general remark, we have also obtained that it is 
actually quite difficult (if not impossible) to fit in a satisfac- 
tory way the simulated rotation curves with Newtonian halo 
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models if the amplitude of the modified potential is set to 
5 = 1, i.e. when the Yukawa term in the point mass case has 
the same weight as the Keplerian one. This is not surprising 
since the boost in the circular velocity in the halo dominated 
regions is so large that can only be reproduced by increas- 
ing the virial mass to unacceptably large values. Although a 
more detailed analysis is needed, this result could be consid- 
ered as an evidence against a too large deviation from the 
Keplerian 1/r scaling. Indeed, should the case S = 1.0 be 
realistic, then the actually observed rotation curves would 
resemble the simulated ones and hence we should have been 
unable to fit them. This is obviously not the case since a 
plethora of successful fits are available in literature. We can 
therefore argue that the case 5 = 1 is not realistic at all or, 
in other words, that a Yukawa - like deviation from the New- 
tonian potential must be only a subdominant correction to 
the 1/r scaling in the inner galaxy regions. 

Although a more detailed analysis is needed, we would 
finally stress that our analysis points towards a new usage of 
the rotation curves. Looking for inconsistencies rather than 
for agreement between these data and Newtonian models 
can indeed tell us not only whether the dark matter particles 
properties should be modified or not, but also whether our 
assumptions on the underlying theory of gravity are correct 
or not. Although it is likely that a definitive answer on this 
question could not be achieved in this way, the analysis of the 
rotation curves data stands out as a new tool to deal with 
modified gravity at scales complementary to those tested 
by cosmological probes. Asking for consistency among the 
results on such different scales could help us to select the 
correct law governing the dominant force of the universe. 
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ABSTRACT 

A viable alternative to the dark energy as a solution of the cosmic speed up 
problem is represented by Extended Theories of Gravity. Should this be indeed the 
case, there will be an impact not only on cosmological scales, but also at any scale, 
from the Solar System to extragalactic ones. In particular, the gravitational potential 
can be different from the Newtonian one commonly adopted when computing the 
circular velocity fitted to spiral galaxies rotation curves. Phenomcnologically modelling 
the modified point mass potential as the sum of a Newtonian and a Yukawa -like 
correction, we simulate observed rotation curves for a spiral galaxy described as the 
sum of an exponential disc and a Navarro - Frenk - White (NFW) dark matter halo. We 
then fit these curves assuming parameterized halo models (either with an inner cusp or 
a core) and using the Newtonian potential to estimate the theoretical rotation curve. 
Such a study allows us to investigate the bias on the disc and halo model parameters 
induced by the systematic error induced by forcing the gravity theory to be Newtonian 
when it is not. As a general result, we find that both the halo scale length and virial 
mass are significantly overestimated, while the dark matter mass fraction within the 
disc optical radius is typically underestimated. Moreover, should the Yukawa scale 
length be smaller than the disc half mass radius, then the logarithmic slope of the 
halo density profile would turn out to be shallower than the NFW one. Finally, cored 
models are able to fit quite well the simulated rotation curves, provided the disc 
mass is biased high in agreement with the results in literature, favoring cored haloes 
and maximal discs. Such results make us argue that the cusp/core controversy could 
actually be the outcome of an incorrect assumption about which theory of gravity 
must actually be used in computing the theoretical circular velocity. 

Key words: dark matter - gravitation - galaxies : kinematic and dynamics 



1 INTRODUCTION 

The picture of a spatially flat universe undergoing acceler- 
ated expansion is the nowadays accepted view of our cosmo. 
According to the successful concordance ACDM model (Car- 
roll et al. 1992; Sahni & Starobinski 2000), there are two 
main ingredients in this scenario, namely dark matter (ac- 
counting for the clustering of the structures we observe) and 
the cosmological constant A (dominating the energy bud- 
get and driving the cosmic speed up). The anisotropy spec- 
trum of cosmic microwave background (de Bernardis et al. 
2000; Komatsu et al. 2010), the galaxy power spectrum with 
the imprinted baryon acoustic oscillations (Eisenstein et al. 



2005; Percival et al. 2010) and the Hubble diagram of Type 
la Supernovae (Kowalski et al. 2008; Hicken et al. 2009) 
represent an incomplete list of the wide amount of data this 
model is able to excellently reproduce in a single scenario. 

From a theoretical point of view, however, the ACDM 
is far to be satisfactory. First, the A term is ~ 120 orders 
of magnitude larger than what expected from quantum field 
theory, while the ratio between its energy density and the 
matter one is coincidentally of order unity just today while 
it should have been much smaller or much larger than 1 
over the rest of the universe history. Motivated by these 
unpleasing shortcomings, a plethora of models giving rise 
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to a varying A - like term has been proposed mainly based 
on a scalar field evolving under the influence of its own self 
interaction potential. Needless to say, the absence of any 
candidate for this scalar field and the full arbitrariness in 
the choice of the potential are serious drawbacks of these 
models which therefore represent only a way to change the 
problems without actually solving them. 

On galactic scales, the A term gives a negligible con- 
tribution to the gravitational potential so that the classical 
Newtonian theory is usually adopted. While the evolution of 
the universe is driven by the cosmological constant, the for- 
mation of structure is mainly determined by the dark matter 
(DM) which provide the potential wells where baryons col- 
lapse to originate the visible component of galaxies. Numer- 
ical simulations allow to follow this process predicting the 
structure of DM haloes. Surprisingly, the theoretical expec- 
tations are not in agreement with observations on galactic 
scales. In particular, the density profile of DM haloes is ex- 
pected to follow a double power - law with p tx x~ a (l+x) 3 ~ a 
with x = r/r s , r s a characteristic length scale and a giving 
the logarithmic slope in the inner regions. Although the de- 
bate on what the precise value of a is remains open, what is 
granted is that a should definitely be positive with a — 1 for 
the most popular NFW model (Navarro et al. 1997). The 
DM haloes are thus described by cusped models or, in other 
words, the logarithmic slope 7 = din p/dlnr never vanishes. 
On the contrary, rotation curves of low surface brightness 
galaxies (LSB) are definitely better fitted by cored models, 
i.e. 7 = in the inner regions, such as the pseudo - isothermal 
sphere, p cx 1/(1 + x 2 ) (Binney & Tremaine 1987) or the 
Burkert model, p cx (1 + a;) _1 (l + x 2 )' 1 (Burkert 1995; 
Salucci & Burkert 2000). As well reviewed in de Blok et al. 
(2010), cored models turn out to be statistically preferred 
over cusped ones from the fit to large samples of LSB rota- 
tion curves. Moreover, for those galaxies where a statistically 
acceptable fit for a cusp model is obtained, it turns out that 
the concentration parameter (defined later) is significantly 
larger than what expected from numerical simulations (see 
de Blok 2010 and refs. therein for further details). 

From the above picture, it is clear that, although ob- 
servationally successful on cosmological scales, the ACDM 
model is far to be free of problems, especially if one also 
remember that there is up to now no laboratory final ev- 
idence for any of the many particles candidate to the role 
of DM. It is therefore worth going beyond the usual view 
and look at a radical revision of the underlying scheme. To 
this end, one can consider cosmic speed up as the first evi- 
dence of a breakdown of General Relativity (GR) as we know 
it. Rather than being due to a new actor on the scene, the 
accelerated expansion can be a consequence of gravity work- 
ing in a different way than GR predicts. This consideration 
has attracted most interest towards braneworld - like models 
such as the five dimensions DGP models (Dvali et al. 2000; 
Lue & Starkman 2003) or fourth order theories of gravity, 
where the GR Einstein - Hilbert Lagrangian is generalized 
by the introduction of a function f(R) of the scalar curva- 
ture (Capozziello 2002; Nojiri & Odintsov 2007; Capozziello 
& Francaviglia 2008; Sotiriou & Faraoni 2010; de Felice & 
Tsujikawa 2010; Capozziello & Faraoni 2010). It is worth 
stressing that, notwithstanding which is the correct modi- 
fied gravity theory, should GR be incorrect on cosmological 
scales, one has to check whether the gravitational potential 



on galactic scale is. Most of modified theories indeed induce 
negligible changes to the gravitational potential on Solar 
System scale in order to pass the classical tests of gravity 
(Will 1993), but this does not prevent to have significant 
deviations from the Newtonian potential on the much larger 
scale of galaxies where no direct experimental test is avail- 
able. 

Although modified gravity theories are investigated at 
cosmological scales as alternatives to dark energy, we are 
here interested in whether they can also impact the estimate 
of dark matter properties on galactic scales. Should indeed 
the gravitational potential be modified, the computation of 
the rotation curve must be done in this modified framework. 
On the contrary, one typically assumes that Newtonian me- 
chanics holds and then constrains the halo model parameters 
by fitting the theoretical rotation curve to the observed one. 
We are here interested in investigating whether this incor- 
rect procedure biases in a significant way the determination 
of the halo parameters and whether such a bias can explain 
the inconsistencies among theoretical expectations and ob- 
servations. In a sense, we are wondering whether modified 
gravity could be a possible way to solve the cusp/core and 
similar problems of the DM scenario. 

The plan of the paper is as follows. In Section 2, we 
present the modified gravitational potential used and de- 
tail the derivation of the rotation curve. It is worth noticing 
that Yukawa-like corrections emerges in any analytical f(R)- 
gravity model, except f(R) = R, where R is the Ricci scalar 
adopted in the Hilbert-Einstein action. Section 3 describes 
how we estimate the bias induced on the halo model pa- 
rameters by fitting data with the Newtonian potential while 
the underlying theory of gravity is non-Newtonian. The re- 
sults of this analysis are discussed in Sections 4 and 5, while 
Section 6 summarizes our conclusions. 



2 YUKAWA-LIKE GRAVITATIONAL 
POTENTIALS 

As it is well known, Newtonian mechanics is the low energy 
limit of GR. Indeed, looking for a stationary and spheri- 
cally symmetric solution of the Einstein equations gives the 
Schwarzschild metric whose tt component gives to the New- 
tonian 1/r gravitational potential in the weak field limit. 
Modifying GR leads to modified field equations hence to a 
different solution and potential in the low energy limit. As 
such, we should first choose a modified theory of gravity to 
finally get the modified potential. On the contrary, we adopt 
here a more general approach postulating that the gravita- 
tional potential generated by a point source m is given by : 

^)=-(i?^h^xp(-i)] a) 

where (S, A) depend on the parameters of the theory. It 
is worth noticing that a Yukawa- like correction has been 
invoked several times in the past. For instance, Sanders 
(1984) showed that the observed flat rotation curves of spi- 
ral galaxies may be well fitted by this model with no need 
for dark haloes provided 8 < and A is adjusted on a 
case -by -case basis. More recently, Yukawa- like corrections 
have been obtained, as a general feature, in the framework 
of /(i?)-theories of gravity (Capozziello et al. 2009a) and 
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successfully applied to clusters of galaxies setting 5 = 1/3 
(Capozziello et al. 2009b). In general, one can relate the 
length scale A to the mass of the effective scalar field intro- 
duced by the Extended Theory of Gravity 1 . The larger is 
the mass, the smaller will be A and the faster will be the 
exponential decay of the correction, i.e., the larger is the 
mass, the quicker is the recovering of the classical dynam- 
ics 2 . Eq.(l) then gives us the opportunity to investigate in 
a simple and unified way the impact of a large class of mod- 
ified gravity theories, among these the Extended Theories, 
since other details do not have any impact on the galactic 
scales we are interested in. 

Eq.(l) is our starting point for the computation of 
the rotation curve of an extended system. To this end, we 
first remember that, in the Newtonian gravity framework, 
the circular velocity in the equatorial plane is given by 
v 2 (R) — Rd&/dR\ z= o, with "!> the total gravitational poten- 
tial. Thanks to the superposition principle and the linearity 
of the point mass potential on the mass m, this latter is com- 
puted by adding the contribution from infinitesimally small 
mass elements and then transforming the sum into an inte- 
gral over the mass distribution. For a spherically symmetric 
body, one can simplify this procedure invoking the Gauss 
theorem to find out that the usual result v a (r) = GM(r)/r 
with M(r) the total mass within r. However, because of the 
Yukawa -like correction, the Gauss theorem does not apply 
anymore and hence we must generalize the derivation of the 
gravitational potential 3 . Alternatively, one can also remem- 
ber that v c (R) — RF{R,z = 0) being F the total gravita- 
tional force and R a radial coordinate. This is the starting 
point adopted in Cardone et al. (2010) where a general ex- 
pression has been derived for the case of a generic potential 
giving rise to a separable force, i.e. : 



F p (»,r) = ^fMf r (v) 



with n = tu/Mq, r\ — r/r s and (Mq,7- s ) the Solar mass 
and a characteristic length of the problem. In our case, it is 
f fi = 1 and : 



Mv) 



1 + -^ 



rj_ \ exp(-? 7 /r? A ) 



(1 + 5)t7 2 



(3) 



with rj\ = A/r s . Using cylindrical coordinates (R,8,z) and 
the corresponding dimensionless variables (77, 8, (with ( = 
z/r 3 ), the total force then reads : 



F(r) = 



GpoTs 



1 Referring to /(i?)-gravity, we prefer to deal with "Extended 
Theories of Gravity" and not modified theories of gravity since 
these theories are nothing else but a straightforward extension of 
GR where the considered action is a generic function of the Ricci 
scalar (Capozziello & Francaviglia 2008; Capozziello & Faraoni 
2010). 

2 Note that the factor 1/(1 + 5) has been explicitly introduced in 
order to recover the correct Newtonian potential for A — > 00. 

3 The non-validity of the Gauss theorem is not a shortcoming 
since, as discussed in Capozziello et al. (2007), physical conserva- 
tion laws are guaranteed by the Bianchi identities that must hold 
in any modified theories of gravity. 



rfdn' / d(' / f r (A)p( v ',e',(')de' (4) 



with p = p/po, po a reference density, and we have defined 



A = [n + V 2 - 2W cos {8 - 8') + (C - C') 2 ] 



/ N 2-|l/2 



{->) 



Since we will be interested in axisymmetric systems, we can 
set p = p{ri,C,). Moreover, the systems of interest here are 
spiral galaxies which we will model as the sum of an infinites- 
imally thin disc and a spherical halo, so that a convenient 
choice for the scaling radius r s will be the disc scale length 
Rd- Under these assumptions, the rotation curve may then 
be evaluated as : 



v 2 c (R) 



with 



1 + 5 



r,'dr,' / p(r/,C')< / MA )d6' (6) 



A = A{8 = C = 0) = [rj 2 + if - 27777' cos 6' + C' 2 ] 



/2]l/2 



(7) 



Inserting Eq.(3) into Eq.(6) gives rise to an integral which 
has to be evaluated numerically even for the spherically sym- 
metric case. It is, however, clear that the total rotation curve 
may be splitted in the sum of the usual Newtonian one and 
a corrective term disappearing for A — » 00, i.e. when the ex- 
tended gravity has no deviations from GR on galactic scales. 



3 ESTIMATING THE BIAS 

Let us assume that a spiral galaxy can be modelled as the 
sum of an infinitesimally thin disc and a spherical halo and 
denote with p the halo model parameters. We can then write 
the total rotation curve as : 



2 c (R,M d , Pl ) = v dN (R,M d )+v 2 hN {R, Pl ) 
+ v 2 dY (R,M d ) + v 2 hY (R, Pi ) 



where M d is the disc mass, the labels d and h denote disc and 
halo related quantities, while TV and Y refer to the Newto- 
nian and Yukawa- like contributions. These latter terms fade 
away for r >> A so that the outer rotation curve is likely the 
same as the Newtonian one. On the other hand, in the in- 
ner region, the two curves may differ more or less depending 
on the value of X/R d . It is worth wondering whether such 
a difference may be compensated by adjusting the model 
parameters. That is to say, we are looking for a new set 
(M' d ,p') such that: 

v 2 N (R, M d , p) = v 2 dN {R, M' d ) + v 2 hN (R, p') . 

Formally, this problem could be solved explicitly writing 
down the above relation for + 1 values of r with Af 
the number of halo parameters and then checking that the 
matching between the two curves is reasonably good (if not 
exact) along the full radial range. Actually, such a procedure 
is far from being ideal since it introduces a dependence of 
the results on the values of r chosen. Moreover, we do not 
need to exactly match the two curves, but only find (M' d , p') 
in such a way that the two curves trace each other within 
the typical observational uncertainties. 
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In order to find (M^,p') taking care of typical observa- 
tions, we therefore adopt the procedure sketched below : 

i. Compute the theoretical circular velocity for input 
model parameters (Md, p) using the exact expression. 

ii. Generate a simulated rotation curve by sampling the 
above v c and adding noise to the extracted points. 

iii. Fit the simulated curve with the same disc + halo 
model, but using the Newtonian theory to compute v c (R). 

iv. Compare the output parameters (M' d ,p') with the 
input ones (Md,p) as function of the normalized scale 
length 77A = X/Rd of the modified gravity theory and other 
quantities of interest. 

In the following, we describe in more details steps i. and ii., 
while Sections 4 and 5 are devoted to discussing the results 
from a large sample of simulated curves. 



3.1 Modeling spiral galaxies 

As yet said above, we will model a spiral galaxy as the sum 
of a thick disc and a spherical halo. For the disc, we adopt 
a double exponential disc so that the density profile reads : 



g(x) = ln(l+x) - x/(l + x) 



(11) 



Pd{R,z) = 



AmRZzd 



exp 



R_ 

Rd 



(8) 



where (Md,Rd,Zd) are the disc mass, lengthscale and 
heightscale. The Newtonian rotation curve cannot be com- 
puted analytically, but, if Zd « Rd (as we will assume), it 
is well approximated by the formula for the infinitesimally 
thin disc (Freeman 1970; Binney & Tremaine 1987) : 



v 2 dN (R) = (2GM d /Rd)y 2 [Io(y)K (y) - Ji(y)ffi(y)] 



(9) 



with y — R/2Rd and I n {y), K„(y) the modified Bessel func- 
tions of order n of the first and second kind, respectively. 

While the observed photometry motivates the use of 
the exponential profile for the disc, the choice of the dark 
halo model is not trivial. Numerical simulations of struc- 
ture formation are typically invoked as a direct evidence 
favouring the use of the NFW density law or its variants. 
However, the NFW model is the outcome of DM only sim- 
ulations performed in a Newtonian framework, while here 
we are working in a modified gravity theory. In principle, 
one should therefore rely on the results of simulations which 
include both the effect of the different potential and the im- 
pact on the evolution of structure due to deviations from 
GR. To this end, one has first to specify which is the mod- 
ified gravity theory one is considering, i.e., explicitly write 
down the gravity Lagrangian. In the case of /(-R)-gravity, 
Schmidt et al. (2009) have shown that, provided f(R) sat- 
isfies the constraints from the Solar System, the halo density 
profiles of DM haloes from numerical simulations is still well 
approximated by the NFW model over the range of masses 
of interest here. Motivated by this result, we therefore as- 
sume that the DM density profile is the NFW one : 

A/lvir 



Ph{r) 



with 



A-KR 3 a g(R vir /R a 



(10) 



In Eq.(10), M v ir and R V i r are the virial mass and radius. 
They are not independent being related by 

1/3 



Rvir — 



/ 3M„i 



\4ivAthPM 



with A t h the overdensity for spherical collapse and pm = 
3H$tl M /&nG the mean matter density today. We follow 
Bryan & Norman (1998) for At;, and set (Q.M,h) = 
(0.28,0.70) in accordance with Komatsu et al. (2010). 

Because of the spherical symmetry, the Newtonian ro- 
tation curve may be easily evaluated as : 

GM h {r) _ GM mr g(r/R s ) 



v h N{r) 



(12) 



r Rvir g(Rvir/Rs) 

While the Newtonian contributions to the rotation curve 
may be computed in the usual way, the Yukawa -like terms 
in the potential give rise to two further terms in v 2 (R) that 
have to be computed numerically. To this end, we must only 
insert Eqs.(8) and (10) into Eq.(6) with / r (Ao) given by 
Eq.(3) subtracting the first term in parentheses. Some care 
must be taken in choosing the reference radius r s . Since 
the data typically probe a limited range in Rd, a natural 
choice is to set r s — Rd. However, the reference density po 
is not the density at r a . It is indeed more convenient to 
set po = pd(Rd,0) for the disc and po ~ Ph{Rs) for the 
halo. With such a choice, for the halo, the dimensionless 
density profile entering Eq.(6) is given by the r- dependent 
part of Eq.(10) provided r/R s is replaced everywhere by 
T)/rj s - Finally, we remember the reader that, when computing 
the total rotation curve, the Newtonian terms, for both the 
disc and the halo, given by Eqs.(9) and (12) respectively, 
must be rescaled by the factor 1/(1 + S) in order to recover 
the classical results in the GR limit (A — ¥ oo). 

3.2 Simulating the rotation curve 

In order to be useful, our approach should rely on simulated 
rotation curves that are as realistic as possible. By this, we 
mean that a.) they must refer to spiral galaxies with reason- 
able values of the model parameters and b.) the sampling 
and the noise should be the same as actual data. Point a.) is 
the easiest to address. As a first step, we randomly generate 
the disc scalelength Rd and the halo virial mass from flat 
distributions over the ranges : 

0.5 < R d /Rd,Mw < 2.0 , 11.5 < log M vir < 13.5 , 

with Rd,MW = 2.55 kpc the disc scalelength for the Milky 
Way (Dehnen & Binney 1998; Cardone & Sereno 2005). 
In order to set the disc mass, we first define the halo mass 
fraction within the optical radius as : 

M h (R opt ) 



}dm 



(13) 



M d + M h (R opt ) 

with R pt = 3.2Rd the optical radius and we have approx- 
imated the disc mass within Rd with the total disc mass. 
We then randomly generate /dm from a flat distribution in 
the range (0.9, 1.1)/dm,/« and /dm, fid = 50% (see, e.g., 
Williams et al. (2010) and refs. therein). The halo scale- 
length R s is computed as R a = Rvir I c where the concen- 
tration c is randomly generated from a Gaussian centred 
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Figure 1. Examples of simulated rotation curves with superimposed theoretical curves. From left to right, model parameters are 
{log M d ,logM vir ,c,f dm Aogrix) = (11.15,12.90,10.24,0.47,0.36), (10.90,11.76, 14.77,0.45,-0.92), (10.04,12.10,13.76,0.54,1.11), while 
the simulation parameters are set as discussed in the text. Note that, depending on how the model parameters are set, it is possible to 
get rotation curves which are flat, decreasing or increasing in the outer region. 



16.7 



M vi 



10 11 h- 1 Mr. 



(14) 



and variance set to 10% of the mean. Note that the above 
relation has been derived by Napolitano et al. (2005) for 
the mass range (0.03, 30) x 10 12 M Q following the method de- 
tailed in Bullock et al. (2001) and updating the cosmological 
model. Finally, we need to set the modified potential param- 
eters (8, A). Following the result obtained for /(i?)-gravity 
(Capozziello et al. 2009b), we first set 5 = 1/3 and run dif- 
ferent simulations randomly generating logr]\ = log (X/Rd) 
from a flat distribution covering the wide range (—2,2). In 
order to explore the impact of S, we also consider the ex- 
tremal case S = 1.0 thus maximizing the contribution of the 
Yukawa terms. 

Having thus generated a realistic galaxy model, we now 
need a rule for sampling it and adding noise in such a way 
that the simulated rotation curve is similar to an observed 
one. A unique choice is not possible since the details of any 
observation depend on the instrumental setup, the observing 
conditions and the galaxy surface brightness. After a visual 
examination of rotation curves samples in literature, we have 
adopted the strategy summarized below. 

(i) Take 2M s im equally spaced points r]i in the range 

r)max)- and replace each of them with rji = £i?7i with 
Ei a randomly generated from the range (0.9, 1.1). 

(ii) For each point in the sample, generate v s i m (j\i) from 
a Gaussian distribution centred on the theoretical value 
and with variance set to (e c /'2)v c (fji). 

(iii) Set the error on the i-ih point as a% = 5ie c v s im(fji) 
with Si randomly chosen in the range (0.9, 1.1). 

(iv) Generate two random numbers (ui,ua) in the range 
(0, 1) and take the point i if ui < Ui- 



(Afsim, Vmin, Vmax,£c) = (30, 0.1, 3.1, 0.25) 

and then 

(A4 ) = (20,3.0,10.0,0.20) . 

We then add the two samples and rescale the errors in such 
a way that the standard \ 2 equals 1 for the input rotation 
curve. Although such values are arbitrary, we have checked 
that the simulated rotation curves have a similar sampling 
and uncertainties of many dataset in literature (see Fig. 1) 
so that we will not try other possible combinations. In or- 
der to quantify the bias on the model parameters, we then 
simulate 1000 rotation curves and fit different models (as 
described in the following) to determine their best fit pa- 
rameters. Note that we will not consider the errors on the 
fit quantities since they strongly depend on the uncertain- 
ties on the data points hence on the choice of the simulation 
parameters (Af B im,,rjmin,rjmax,£ c )- Since we are interested 
in a statistical analysis of the full sample rather than on a 
detailed study of a particular case, we prefer to limit our 
attention to the best fit parameters only thus avoiding to 
deal with the details of the simulation procedure. 



4 BIAS ON HALO PARAMETERS 

The sample of simulated rotation curves constructed by the 
procedure detailed above is the starting point of our anal- 
ysis. Indeed, we can now fit each curve with a given (New- 
tonian) model and compare the output best fit parameters 
with the input ones. There is, however, a preliminary caveat 
to be addressed. When fitting a whatever model to a given 
dataset, one has to put down an objective criterium to deem 
the model as reliable or not. It is then common practice to 
compute the standard \ 2 , i.e. 



2 

X = 



Af r 

E 



Vc(Ri) - V c>t h(Ri) 



Typical rotation curves are well sampled up to the optical 
radius R op t = 3.2Rd, while the sampling gets worse at larger 
radii. On the contrary, the percentage errors are of the same 
order in both regions, although it is possible that the inner- 
most points have larger uncertainties because of deviations 
from ordered motions. For given input model parameters, 
we therefore generate two samples of points by first setting 



and then evaluate the quality of the fit considering the value 
of x 2 = X 2 /d.o.f. with d.o.f. = AT — Af p the number of de- 
grees of freedom and Af (Af p ) the number of data points (pa- 
rameters). Formally, one can then conclude that the model 
is a good description of the data if \ 2 ^ Xth with the thresh- 
old value depending on Af. Since, for our simulated curves 
Af ~ 50, Xth ~ 1-2 which is highly demanding selection 
criterium. Actually, this formal rule rigorously applies only 



© 0000 RAS, MNRAS 000, 000-000 



6 V.F. Cardone & S. Capozziello 



if the data points are uncorrelated and the measurement 
uncertainties are Gaussian distributed. Both these assump- 
tions break down for our simulated curves since we have 
rescaled the uncertainties in such a way that x 2 = 1 for 
the input model. By a visual examination of many fitted 
curves, we have checked that a reasonably good fit (with 
no systematic deviations from the outer or inner data and 
rms percentage deviations smaller than 10%) is obtained up 
to x 2 ~ 2.5 which we choose as our cut. We will therefore 
consider a model as successfully fitting the data if x 2 ^ 2.5 
and log M v i r ^ 14 with this latter criterium introduced to 
avoid unphysical solutions. However, we will also discuss the 
trends for the quantities of interest for a subsample obtained 
by using a stringent criterium, i.e. x 2 ^ 1.5 in order to ex- 
plore the impact of the threshold x 2 value. 

4.1 Navarro-Frenk- White models 

In realistic situations, one has a set of (R, v c ) data and a 
measurement of the galaxy surface brightness in a given 
band. It is then common to set the disc scalelength Rd to the 
value inferred from photometry, while the disc mass can be 
inferred from the total luminosity provided an estimate of 
the stellar M/L ratio is somewhat available (e.g., from stel- 
lar population synthesis models fitted to the galaxy colours). 
As a first step, we therefore assume that both (Rd, Md) are 
known and fit the simulated data for each rotation curve to 
determine the NFW model parameters 4 only. 

4-1.1 5 — 1/3 simulated samples 

Let us first consider the results obtained setting 5 = 1/3 
when simulating the rotation curves. As yet said above, we 
are here trying to fit modified gravity rotation curves us- 
ing theoretical models computed in a Newtonian framework 
so that the first point to address is whether this is indeed 
possible. To this end, it is sufficient to check what is the per- 
centage of rotation curves passing the two selection criteria 
quoted before. We indeed find that 90% of the simulated 
rotation curves are fitted with x 2 ^ 2.5 and logM„ ir < 14 
(which we will refer to as the well fitted sample, hereafter 
WS), while this fraction decreases to 77% if x 2 ^ 1-5 is de- 
manded (defining what we will hereafter refer to as the best 
fitted sample, BS). Averaging over the sample values gives : 

(X 2 > = 1-27 ± 0.24 (1.19 ± 0.13) for WS (BS) 

while for the rms of the percentage deviations Av c /v c = 
(v° bs - vl u )/vf s we get : 

(rms(Av c /v c )) = 6.4% ± 0.8% (6.3% ± 0.7%) for WS (BS) . 

Motivated by the high fraction of successful fits and the 
low values of x 2 and rms(Av c /v c ), we can therefore safely 

4 To this end, we use a straightforward Monte Carlo Markov 
Chain algorithm to minimize the x 2 with respect to the param- 
eters (log r/s, log V v i r ) and then infer the values of (c, M v i r ) and 
the dark matter mass fraction foM- This approach is more stable 
than fitting directly for (c, M v { r ) and allows to correctly recover 
the input model parameters when we fit rotation curves simulated 
in the Newtonian framework. 



conclude that it is possible to fit a modified gravity rota- 
tion curve with a NFW Newtonian one obtaining both an 
acceptable x 2 value and reasonable model parameters. 

Having checked that the model reasonably well fit the 
data, we can rely on the best fit model parameters and com- 
pare them with the input ones to estimate what is the bias 
induced by an incorrect assumption of the gravity law. To 
quantify this bias, we define 1Z(x) = Xf it /x 3 i m , i.e. the ratio 
between the best fit and the input values of a given quan- 
tity x. Needless to say, should IZ(x) be on average close to 
1, we can conclude that assuming Newtonian gravity to fit 
modified gravity rotation curves does not induce a signifi- 
cant bias on the estimate of x. The results, summarized in 
Table 1 for both the WS and BS samples, show that such 
a bias is indeed quite important and only mildly depend on 
the selection criteria adopted (so that we will hereafter refer 
to the WS results only unless otherwise stated). 

Table 1 gives some statistics on the distribution of IZ(x) 
for the halo parameters (r) s ,M v i r ). However, these values 
could be misleading since the histograms for lZ(rj„) and 
lZ(M V i r ) are strongly asymmetric with long tails towards 
the right. In other words, TZ(r/ a ) and lZ(M v i r ) could be much 
larger than their median value with TZ(ri a ) being as high 
as 12 and values of lZ(M v i r ) as large as 15 so that both 
the halo scalelength and virial mass can be grossly over- 

1/3 

estimated. Since c oc R vir oc M v ' ir , one could naively ex- 
pect that the concentration is overestimated too. On the 
contrary, the distribution of the IZ(c) values is almost sym- 
metric around its mean clearly disfavouring IZ(c) > 1, i.e. 
the concentration is underestimated. Actually, such a re- 
sult can be understood remembering that c oc R^ 1 so that 
lZ(c) oc (M vir ) /lZ(r/ a ) finally leading to values smaller 
than unity. The emerging picture is therefore that of a halo 
having a larger mass than the input one, but also a larger 
scalelength. This can be qualitatively explained noting that 
the Yukawa- terms in the rotation curve increases the net 
circular velocity both in the inner and outer regions probed 
by the data. In order to adjust the fit in the halo dominated 
regions, one has to increase M v i r , but then R a has to be 
increased to in order to lower the density (and hence the 
contribution to the rotation curve) in the inner regions not 
to overcome the observed circular velocity. As a result, the 
concentration is underestimated and the same takes place 
for the dark matter mass fraction within the optical radius 
since the halo mass is now pushed outside the optical ra- 
dius thus explaining why TZ(foM) is smaller than unity. As 
Fig. 2 shows, lZ(x) is correlated with log r]\ with the sign of 
the correlation depending on 1Z(x) being typically larger or 
smaller than 1. Not surprisingly, the larger is the Yukawa 
scalelength A, the closer is IZ(x) to 1, i.e., the smaller is the 
bias induced by the assumption of Newtonian gravity. How- 
ever, there is a large scatter in these correlations since the 
relative importance of the correction also depends on the 
input halo parameters. Indeed, when R a ~ A, the halo con- 
tribute to the modified rotation curve is maximized so that 
the Yukawa corrections are more important and the bias of 
the fitted parameters is larger. As a consequence, one can 
therefore expect a correlation between 7Z(x) and log j) s which 
is what we indeed find looking at the results in Table 1. 
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Figure 2. 7Z(x) vs log^ with x = r] s (upper left), M„; r (upper right), c (lower left), Jum (lower right) for the NFW model fit to the 
simulated rotation curves with 8 = 1/3. Note that only 10% of the WS sample is plotted to not clutter the figure. 



Table 1. Bias 1Z(x) on the NFW model parameters assuming the disc mass is known and 8 = 1/3. Columns are as follows : 1. parameter 
id; 2., 3., 4. mean ± standard deviation, median and rms values, 5., 6., 7., 8., 9. Spearman rank correlation coefficients between IZ(x) 
and input logjjs, logM„i r , c, /dm, logrjx- Upper half of the table is for the WS sample, while lower half for the BS one. 



Id 


{11) 


Turned 


^■rms 


C(logr, 3 ,n) 


C(\ogM vir ,Tl) 


C(c, Tl) 




c{\ ognx ,n) 


Is 


4.0 ±2.8 


2.9 


4.8 


0.31 


0.18 


-0.17 


-0.09 


-0.47 




3.9 ±3.8 


2.4 


5.5 


0.48 


0.31 


-0.27 


-0.05 


-0.32 


c 


0.45 ±0.15 


0.45 


0.48 


-0.13 


-0.05 


0.07 


0.13 


0.55 


Sou 


0.70 ±0.06 


0.70 


0.71 


0.36 


0.30 


-0.20 


0.39 


0.70 


Us 


3.8 ±2.8 


2.8 


4.7 


0.32 


0.24 


-0.23 


-0.07 


-0.42 




3.6 ±3.6 


2.2 


5.1 


0.50 


0.36 


-0.34 


-0.02 


-0.26 


C 


0.47 ±0.15 


0.47 


0.49 


-0.14 


-0.10 


0.13 


0.12 


0.53 


fnu 


0.70 ±0.06 


0.70 


0.71 


0.39 


0.31 


-0.18 


0.40 


0.71 



4-1.2 5 = 1.0 simulated samples 

Let us now consider the results for the rotation curves sam- 
ple simulated under the extreme assumption 5 = 1.0, i.e. 
when the contribution of the Yukawa term to the point mass 
potential (1) is the same order of magnitude as the Newto- 
nian one. Note that, different from the previous case, S = 1.0 
is, as far as we know it, an unmotivated assumption, but we 
consider it to investigate the impact of 5 on the results. 

The first striking result is that the fraction of well fit- 
ted rotation curves dramatically drops to 22% applying the 
WS selection criteria (x 2 ^ 2.5 and logM„i r ^ 14), while no 
curves pass the cut x 2 ^ 1-5 thus showing that it is likely 
not possible to reproduce the simulated curves with NFW 



model if Newtonian gravity is assumed. Such a result can 
be understood noting that, for S = 1.0, the modified grav- 
ity rotation curve is much larger than the Newtonian one 
so that one has to greatly increase the halo virial mass to 
fit the outer rotation curve. In order to compensate the cor- 
responding increase in the inner region, one should set rj B 
as large as possible, but this also tend to lower the circu- 
lar velocity for R >> R 3 . Finding a compromise between 
these two opposite trends is quite difficult thus leading to 
unacceptably large x 2 values. 

It is worth investigating what is the bias induced on 
the halo model parameters for the successfully fitted cases. 
The trends of the bias parameters IZ(x) are the same as 
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Figure 3. Histogram of "R-(a) values and its dependence on \ogn\ for the gNFW model (with disc mass set to the input value) fit to 
the simulated rotation curves with <5 = 1/3. Note that only 10% of the WS sample is plotted to not clutter the figure. 



in the previous case, but the typical values are now much 
larger. In particular, 1Z{ris) can be as high as 40 thus leading 
to grossly underestimates of both the concentration c and 
the dark matter mass fraction /dm- Not surprisingly, the 
only way to reduce such large biases is to increase logr/x 
to values larger than our upper limit log rj\ — 2. However, 
in this case, the modified gravity potential differs from the 
Newtonian one only on group and cluster scales where the 
rotation curve are no more measured so that our analysis 
becomes meaningless. 



4.2 Generalized Navarro-Frenk- White models 

Up to now, we have considered the NFW profile for the 
halo model, but such a choice does not allow us to inves- 
tigate whether the mismatch between modified and Newto- 
nian gravity can have an impact on the inner logarithmic 
slope. In order to explore this issue, we therefore consider 
the generalized NFW model (gNFW, Jing & Suto 2000) : 



Rvir 
~R~> 



Rvir 



(2 - a)R s 



(17) 



Ph(r) = 
with 

h(x) = 



±nRlh{R vir /R s ) \R 



1 + 



R a 



-(3-c0 



(15) 



(16) 



The parameter a depends on the behaviour of the rotation 
curve in the inner regions where the disc contribute can be 
larger than the halo one. It is therefore worth investigating 
which is the impact of the uncertainties on Md in order to 
see how the bias on a depend on it. 



4-2.1 gNFW models with known disc mass 

As a first step, we assume that Md is set to the input value 
and only look at the bias on the parameters (a, c, M v i r ). As 
a preliminary remark, it is mandatory explaining how the 
concentration is defined for gNFW models. To this end, one 
may simply note that, for the NFW density profile, R 3 is the 
radius at which the logarithmic slope equals the isothermal 
value 7 = —2. Therefore, the concentration for the gNFW 
model may be easily defined as : 



with j{R- 2 ) = -2. Note that, for a = 1, the gNFW model 
reduces to the NFW one and Eq.(17) gives the standard 
definition for the concentration of NFW haloes. 

Applying the selection criteria used above to the sample 
with S = 1/3, we find that the WS sample contains now 90% 
of the simulated curves, while this fraction lowers to 82% for 
the BS sample. Averaging over the samples gives : 

(x) = 1-28 ± 0.19 (1.23 ± 0.12) for WS (BS) , 

(rms(Av c /v c )) = 6.4% ± 0.8% (6.3% ± 0.8%) for WS (BS) . 

These values are almost identical to those obtained for the 
fits of the NFW model described before and allows us to 
confidently conclude that it is possible to reproduce the sim- 
ulated rotation curves with a gNFW model and Newtonian 
gravity. In a sense, this is not surprising since the gNFW 
model has one more parameter than the NFW one so that 
there is a much larger freedom to adjust the shape of the 
theoretical rotation curve to fit the observed one. As a re- 
sult, the percentage of well fitted rotation curves slightly 
increases, but the quality estimator (x 2 and the rms value 
of Av c /v c ) stay almost unchanged. 

The bias induced on the model parameters can be still 
estimated considering the quantity IZ(x) defined before and 
summarized in Table 2 for the present case 5 . Comparing to 
the results in Table 1, it is apparent that the bias induced 
on the halo parameters in common, i.e. (rj 3 , c, M V i r , /dm), 
is almost the same both for the WS and BS samples. The 
same qualitative discussion can be repeated here so that we 
will not consider anymore this issue. It is, on the contrary, 
more interesting to look at the bias on a which asks for some 
caution. As can be seen from the left panel in Fig. 3, the dis- 
tribution of a values is quite large and asymmetric, while the 
right panel is a clear evidence that its mean value strongly 
depend on log n\. Indeed, should we have only fitted models 
with A ^ Rd, IZ(a) would have been smaller than 1, i.e. the 
inner slope would have been underestimated and shallower 



5 Note that, for all the simulated curves, it is a = 1 since we 
have assumed a NFW model in the simulation process. Therefore, 
Tl(a) = a for the gNFW models. 
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Table 2. As in Table 1 for the gNFW with fixed disc mass case. 



Id 


(R) 






C(\ogri s ,n) 


c(iogM vir ,n) 


C{c,K) 


C(/ DM ,K) 


C{\ogri x ,K) 


a 


0.92 ±0.22 


0.95 


0.95 


-0.56 


-0.58 


0.37 


0.09 


0.32 


Is 


2.8 ±2.6 


2.0 


3.8 


-0.11 


-0.06 


0.0 


-0.05 


-0.06 




2.4 ± 2.6 


1.4 


3.5 


0.43 


0.39 


-0.39 


0.03 


-0.18 


c 


0.56 ±0.18 


0.53 


0.59 


0.15 


0.08 


-0.03 


-0.03 


0.16 


Idm 


0.70 ±0.06 


0.71 


0.71 


-0.06 


-0.18 


0.04 


0.40 


0.70 


a 


0.93 ±0.22 


0.96 


0.95 


-0.56 


-0.50 


0.39 


0.14 


0.31 


Us 


2.8 ±2.5 


2.1 


3.7 


-0.08 


-0.07 


0.0 


0.11 


-0.07 




2.4 ± 2.6 


1.4 


2.7 


0.42 


0.36 


-0.37 


0.05 


-0.17 


C 


0.55 ±0.18 


0.53 


0.58 


0.12 


0.08 


-0.03 


-0.09 


0.18 


Jdm 


0.70 ±0.06 


0.70 


0.70 


-0.05 


0.18 


0.03 


0.42 


0.69 



0.4 0.6 0.8 1.0 

K(«) 



Figure 4. Histograms of lZ(x) with x = Md (left) and x = a (right) for the gNFW model with free disc mass fits and <5 = 1/3. 



models be preferred. Should the modified gravity scalelength 
be smaller than the typical disc one, fitting rotation curves 
assuming the validity of Newtonian gravity would have us 
led incorrectly to believe that the inner slope of the density 
profile is much smaller than the NFW one. Such a result, 
therefore, suggests that the cusp/core controversy is just a 
fake problem originated by an incorrect assumption of what 
the underlying gravity theory actually is. 

As a final remark, it is worth noting that a second dif- 
ference among the results in Tables 1 and 2 is represented 
by the markedly different values of the Spearman correla- 
tion coefficients. Actually, this is a consequence of having 
added one more parameter which introduces a degeneracy 
between the fitted quantities hence partially washing out the 
correlations of IZ(x) with the input halo parameters and the 
modified gravity scalelength. In particular, a and log r/ s are 
clearly correlated since both set the shape of the inner rota- 
tion curve. Although with a large scatter, lZ(a) is typically 
smaller (larger) than 1 for log?7A < (> 0) so that TZ(ct) 
actually correlates with logrjx, but not linearly thus giving 
rise to a small Spearman correlation coefficient (which looks 
for linear correlations). Since rj s has to be adjusted accord- 
ing to the changes in a, 7Z(ri 3 ) has to follow TZ(ct) thus losing 
its correlation with log 77a. 

Such results are strengthened if we consider the simu- 
lated curves for S = 1.0, but now relying on a much smaller 
statistics. Indeed, the WS sample is now made out of only 
20% of the full set of simulated curves, while only 2% of 
them are included in the BS sample. These fractions are 



larger than in the NFW fits considered before, but are still so 
small to make us conclude that such extreme modified grav- 
ity curves can not be reproduced in a Newtonian framework 
using the gNFW model. For the few surviving curves, we 
can repeat the same discussion above with the only caveat 
that now the IZ(x) values deviate more from 1 than in the 
5=1/3 case. Moreover, the slope a can be underestimated 
also when logTiA > depending on the input halo parame- 
ters. This is a still further evidence in favour of the suggested 
interpretation of the cusp/core controversy as the outcome 
of a systematic error in the assumed gravity theory. 

4-2.2 gNFW models with free disc mass 

Up to now, we have set the disc mass to the input value im- 
plicitly assuming that one is able to perfectly estimate this 
quantity from the galaxy colors and a stellar population syn- 
thesis code. Actually, this is not always the case and, on the 
contrary, it is usual to include Md in the set of parameters 
to be determined by fitting the rotation curve data. We have 
therefore repeated the above analysis for the gNFW model 
leaving Md free to be adjusted by the fit. 

Considering there is one more fit parameter, it is not 
surprising that the fraction of successfully fitted rotation 
curves (for 5 = 1/3) increases to 92% (85&) for the WS (BS) 
samples and also the quality of the fit is improved being : 

(X 2 ) = 1.19 ± 0.20 (1.14 ± 0.13) for WS (BS) , 

(rms(Av c /v c )) = 6.1% ± 0.6% (6.0% ± 0.8%) for WS (BS) . 
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Table 3. As in Table 2 for the gNFW with free disc mass case. 



1,1 


{'<■/ 


7? 

'Mnefl 


1? 




<s(iogr) s , i-C) 


r*n ntT A/T - "7? \ 




H/DMi ><■> 




M d 


0.79 ±0.19 


0.80 


0.82 


0.0 


0.14 


0.12 


-0.07 


-0.02 


0.09 


a 


1.00 ±0.20 


1.02 


1.01 


-0.10 


-0.01 


-0.08 


-0.01 


0.01 


0.33 




1.6 ± 1.5 


1.2 


1.6 


-0.01 


0.21 


0.18 


-0.15 


0.03 


0.0 




1.7 ±3.5 


0.8 


3.8 


0.06 


0.22 


0.23 


-0.18 


0.02 


-0.09 


C 


0.87 ±0.37 


0.82 


0.95 


0.0 


-0.22 


-0.21 


0.16 


-0.03 


0.07 


Idm 


0.97 ±0.24 


0.96 


1.00 


-0.04 


-0.12 


-0.12 


0.07 


0.05 


0.10 


M d 


0.78 ± 0.20 


0.78 


0.80 


0.03 


0.11 


0.13 


-0.09 


-0.04 


0.08 


a 


1.00 ±0.20 


1.02 


1.02 


-0.12 


0.03 


-0.06 


-0.02 


0.02 


0.35 


Is 


1.6 ± 1.6 


1.2 


2.2 


0.00 


0.20 


0.18 


-0.15 


0.01 


-0.01 




1.6 ±3.2 


0.8 


3.5 


0.07 


0.21 


0.23 


-0.18 


0.0 


-0.12 


C 


0.89 ± 0.38 


0.84 


0.97 


-0.02 


-0.22 


-0.22 


0.17 


-0.01 


0.09 


Sdm 


0.99 ± 0.24 


0.97 


1.01 


-0.06 


-0.09 


-0.12 


0.08 


0.07 


0.10 



Table 3 summarizes the results on the bias estimator 7Z(x) 
for the different parameters involved which allow us to draw 
some interesting considerations. First, we note that the halo 
model parameters are now better recovered than in previous 
cases. For instance, although the corresponding distributions 
have long tails up to 1Z{x) ~ 10 — 15, most of the values of 
7Z(r] s ) and 7i(M V i r ) are smaller than 2 leading to a modal 
value (i.e., the most peak of the histogram) close to 1, i.e. 
(r)s,M v ir) are not biased anymore for most of the WS and 
BS sample curves. A similar discussion also apply to the con- 
centration c and the dark matter mass fraction Jdm whose 
histograms are quite symmetric and centred close to the no 
bias value, 1Z(x) = 1. Concerning a, one must still keep in 
mind that its mean value depend on the cut on log rj\ typi- 
cally being 7Z(a) < 1 (i.e., the inner slope is shallower than 
the input one) when A < Rd- Note, however, that this effect 
is now less pronounced than before thus explaining why the 
distribution of 71(a) values is peaked close to 1 with a not 
too large width. On the contrary, the disc mass M d turns 
out to be biased low with Md being smaller than the input 
value for most of the cases. Such an unexpected result can be 
qualitatively understood considering that the Yukawa terms 
make the modified gravity rotation curve larger than the 
Newtonian one in the inner and outer regions. As in the 
case of the NFW fit, one must increase the halo virial mass 
to recover this additional contribution in the halo dominated 
regions which increases also the halo inner circular velocity. 
In the NFW case, one has then to increase -q s to prevent 
overestimating the simulated curve, while the same effect is 
now accomplished by lowering the disc mass. This explains 
why the halo parameters turn out to be less biased than in 
the NFW fit case. It is, however, worth stressing that, as 
shown in the left panel of Fig. 4, the 7t(Md) distribution is 
actually bimodal. As a consequence, when the disc mass is 
not biased low, one has again to increase rj 3 thus explaining 
the long tails towards 7t(r] a ) » 1 values. Alternatively, one 
can leave rj s unchanged, but make the density profile shal- 
lower hence giving rise to the tail towards small 71(a) < 1 
in the right panel of Fig. 4. Model degeneracies also explain 
why the correlation coefficients among 7Z(x) and the input 
model parameters are typically quite small. We, however, 
caution the reader that a small value of Spearman correla- 
tion does not necessarily imply that 7Z(x) doe not depend 
on the corresponding quantity (see, e.g., the TZ(a) depen- 



dence on log?7A for the gNFW fits with fixed disc mass). 
We have therefore checked whether this is the case or not 
by directly looking at the 7Z(x) vs pi plots (with pi the i- 
th input parameter) finding that the scatter of the points in 
each plot clearly indicates a complicated dependence a weak 
dependence on all the input parameters but log?7A. 

Finally, we have repeated the above analysis setting 
5 = 1.0 to generate the simulated rotation curves. Differ- 
ently from the previous cases considered up to now, we now 
find that it is possible to fit the gNFW model with free disc 
mass to the simulated rotation curves. In a sense, the addi- 
tional freedom we have now to adjust the disc mass makes 
it possible to recover the inner rotation curve better than in 
the previous models thus leading to a successful fit for 97% 
(85%) of the galaxies using the selection criteria for sample 
WS (BS). The quality of the fit is also remarkably good : 

<X 2 ) = 1.32 ± 0.30 (1.18 ± 0.15) for WS (BS) , 

(rms(Av c /v c )) = 6.5% ± 1.0% (6.2% ± 0.7%) for WS (BS) . 

Concerning the values of 7Z(x) quantities, they are almost 
the same as those in Table 3 for the parameters (rj s ,c, /dm) 
which are therefore almost unbiased. On the contrary, for 
the other parameters (and the WS sample), we get : 

0.54 ± 0.20 for x = M d 

(K(x))={ 1.04 ±0.28 for x = a 

0.85 ± 1.18 for x = M vir 
0.58 for x = M d 



[R.(x)] rr , 



1.08 for x = a 



1.5 



for X = M v i r 



As it is clear from the large standard deviations and rms 
values, the distribution of 7Z(M d ) and TZ(M v i r ) are strongly 
asymmetric and, as a result, the one for a presents a non 
negligible amount of points with TZ(a) > 1. This result can 
be explained noting that setting S = 1.0 maximizes the con- 
tribution of the Yukawa terms so that we now have to mimic 
a similar effect by adjusting the model parameters in the 
Newtonian fitting. One possibility is to leave a unchanged 
with respect to the input value (a = 1.0), but increasing 



© 0000 RAS, MNRAS 000, 000-000 



Galaxy haloes and Yukawa-like gravitational potentials 11 



Table 4. As in Table 1 for the PI model with free disc mass case. 



Ill 


(K) 






C (log M d ,K) 


C(logr, s ,ft) 


C{logM vir ,K) 


C(c, 1Z) 


CUdm,K) 


C(log Vx ,K) 


M d 


1.20 ±0.10 


1.20 


1.20 


-0.12 


-0.26 


-0.21 


0.09 


0.55 


0.48 




5.2 ±7.0 


3.8 


8.7 


0.02 


0.12 


0.09 


-0.03 


0.0 


0.06 


Som 


0.37 ±0.08 


0.38 


0.38 


0.08 


0.49 


0.38 


-0.30 


-0.06 


-0.22 


M d 


1.19 ±0.10 


1.19 


1.19 


-0.13 


-0.27 


-0.19 


0.15 


0.57 


0.44 




5.4 ±8.9 


3.7 


10.0 


0.0 


0.06 


0.04 


0.0 


-0.03 


0.09 


Sum 


0.38 ± 0.08 


0.39 


0.39 


0.07 


0.56 


0.42 


-0.34 


-0.10 


-0.27 



Table 5. As inTable 1 for the Burkcrt model with free disc mass case. 



Id 


(K) 


Turned 




C(logA/ d ,ft) 


C(logr, s ,^) 


c{\o g M vir ,n) 


C(c, TV) 


c(f DM ,n) 


C(\og Vx ,1Z) 


M d 


1.14 ±0.11 


1.13 


1.15 


0.0 


-0.09 


-0.03 


-0.02 


0.32 


0.30 




1.1 ± 1.1 


0.7 


0.8 


0.06 


0.13 


0.16 


-0.15 


0.09 


0.13 




0.46 ±0.11 


0.44 


0.48 


-0.06 


-0.03 


-0.07 


0.11 


0.05 


-0.08 


M d 


1.11 ±0.10 


1.10 


1.12 


-0.11 


0.03 


0.05 


0.05 


0.33 


0.49 




1.0 ± 1.1 


0.6 


1.5 


0.05 


0.12 


0.14 


-0.10 


-0.10 


0.21 




0.48 ±0.11 


0.47 


0.50 


-0.01 


-0.05 


-0.07 


-0.01 


0.11 


-0.22 



the halo mass by a significant factor thus leading to large 
values of 1Z(M V ir)- But, in this case, one has also to decrease 
the disc contribution to the inner rotation curve to compen- 
sate for the additional dark matter. On the contrary, one 
can leave almost unchanged M V i r , but redistribute the dark 
mass pushing it towards the inner region. To this aim, one 
should make the profile steeper, i.e. increase a, while M d 
must be smaller to not overcome the inner circular veloc- 
ity. This strategy will lead to TZ(a) > 1 and lZ(M V i r ~ 1, 
while again it is lZ(M d ) < 1. The final histograms of 1Z(x) 
for (M d , a, M v i r ) is then the outcome of a mixture of these 
two possible strategies whose relative contribution depends 
on the details of the individual rotation curves. 



5 BIAS ON THE HALO DENSITY PROFILE 

The use of the gNFW model to fit the simulated rotation 
curve may also be read as a first step towards completely 
mismatching the halo density profile. In a sense, the gNFW 
model departs from the input NFW one only because of 
the possibly different inner logarithmic slope. It is, however, 
common in data analysis to use also models that differs from 
the NFW one both in the inner and outer regions. Moreover, 
the gNFW is still a cuspy model, while the cusp/core contro- 
versy comes from the observation that cored models better 
fit the observed rotation curves. To this end, we now drop 
the implicit assumption that the trial model tracks or gener- 
alizes the NFW profile and investigate whether completely 
different models in the Newtonian framework are in accor- 
dance with our modified gravity simulated rotation curves. 

5.1 The pseudo - isothermal sphere 

A classical example of a cored model is represented by the 
pseudo - isothermal (hereafter, PI) model whose density pro- 
file reads : 



Ph{r) = ^mK-(R-jR s [ 1 + R- S ) (18) 

with 

h(x) = x — arctana; . (19) 

As it is clear, the PI density law approaches a constant value 
for r << R s , while falls off as r~ 2 in the outer regions. As 
such, the PI model is radically different from the NFW and 
gNFW ones so that it is interesting to see whether can fit or 
not the data. Following common practice, we leave the disc 
mass as an unknown and estimate (M d ,R s , M v i r ). 

Applying the corresponding selection cuts, we find that 
78% of the simulated galaxies fall into the WS sample, while 
the BS one contains 43% of the full sample. Such high frac- 
tion and the low values of both the reduced x 2 

(X 2 > = 1.49 ± 0.35 (1.24 ± 0.15) for WS (BS) , 

and of the rms of percentage residuals 

(rms(Av c /v c )) = 7.0% ± 0.9% (6.9% ± 0.9%) for WS (BS) 

make us confident that the Newtonian circular velocity of 
the PI model can indeed mimic the rotation curve predicted 
by modified gravity. This is in agreement with what is indeed 
found in the literature where galaxies rotation curves data 
are typically well fitted by PI haloes. 

Having used now a different halo model to fit the data, 
it is not possible to straightforwardly compare the output 
parameters with the input ones. First, although we use the 
same symbol, the scalelength R s of the PI model is not de- 
fined by the condition "/{Rs) = — 2 as for the NFW model, 
but rather represents the core radius. As such, it is meaning- 
less to compare the input rj s with the output R s /R d since 
they refer to a different characteristic of the density profile. 
Moreover, a concentration for the PI model is not defined 
so that we can not compare with the input one. On the con- 
trary, the total disc mass M d , the virial mass M v i r and the 
dark matter mass fraction /dm refer to global properties 



© 0000 RAS, MNRAS 000, 000-000 



12 V.F. Cardone & S. Capozziello 



of the disc + halo model so that it makes sense to compare 
them with the input ones. Table 4 then gives the bias on the 
global parameters (Md, M v i r , }dm) both the WS and BS 
samples (having set 5 = 1/3). It is clear that all these three 
quantities are severely biased with (Md, M v i r ) being overes- 
timated and /dm grossly underestimated. Such a behaviour 
can be explained noting that, once the core radius has been 
set, the only way to increase the Newtonian rotation curve 
is to increase the global masses (Md, M V i r ) thus leading to 
asymmetric distributions both peaked in lZ(x) > 1. It is 
worth noting that, if one ignores that the underlying grav- 
ity theory is not Newtonian, a large disc mass will be in- 
terpreted as the presence of a maximal disc which is in- 
deed the case in many spiral galaxies (Palunas & Williams 
2000). Being both M V i r biased high, one could find the re- 
sult IZ(fnM) « 1 an unexpected nonsense. Actually, one 
has first to note that also the disc mass is overestimated 
which automatically reduces the dark matter mass fraction. 
Moreover, for the PI model, the mass profile increases al- 
most linearly, while, for the input NFW profile, Mh(r) grows 
approximately in a logarithmic way. Being the input R 3 al- 
most the same as the output core radius, we therefore get 
MNFw{Ro P t) > Mpi(R pt) even if the virial mass of the PI 
model is larger than the NFW input one. As a consequence, 
fr>M turns out to be underestimated as we indeed find. 

As a final test, let us consider what happens when 5 = 
1.0 is used in simulating the rotation curves. It turns out that 
the percentage of curves entering the WS sample reduces to 
67%, while only a modest 12% enter the BS sample. This 
sudden drops is actually due to the cut on the output virial 
mass rather than on the reduced \ 2 ■ As such, we believe 
that such a case can not be mimicked by a PI model under 
the Newtonian gravity assumption and not discuss anymore 
the bias on the output parameters. 

5.2 The Burkert model 

Another cored model but with a different outer profile is the 
Burkert model whose density profile read (Burkert 1995) : 

Ph(r) = 4,i8Ml r W + + 2 (20) 

with 

h B {x) = ln(l + x) -arctana; + (l/2)ln(l + a; 2 ) . (21) 

Note that the Burkert model presents an inner core with ra- 
dius R s , but asymptotically drops off as r -3 (hence having 
a finite total mass). As such, it represents a sort of compro- 
mise between the cored PI model and the cusped NFW one. 
Again, we will leave the disc mass free and determine the 
three parameters (Md,ri a , M v i r ) from the fit. 

Setting S — 1/3, we find that 88% of the galaxies fall 
into the WS sample, while this fraction drops to 37% for 
the BS one. The quality of the fit may be judged from the 
following figures of merit : 

(x) = 1-63 ± 0.32 (1.31 ± 0.15) for WS (BS) , 

(rms(Av c /v c )) = 7.3% ± 1.0% (6.4% ± 0.7%) for WS (BS) , 

comparable to the values obtained for the PI fits. It is worth 
noting that the fraction of galaxies in the WS sample is 
larger than for the PI model as a result of the different mass 



profile. Indeed, for the PI model, asymptotically M(r) oc r 
so that v c (r) oc M(r)/r ~ const, while for the Burkert model 
we have a finite total mass and hence a Keplerian fall off. As 
Fig. 1 shows, our simulated rotation curves sample is made 
out of galaxies which can be flat, decreasing or increasing in 
the outer regions depending on the input model parameters. 
It is, of course, quite difficult to fit decreasing curves with the 
PI model which does not predict such a behaviour, although 
the limited radial range probed and the uncertainties allow 
to get a reasonable match in some cases. Models with ph oc 

for r >> R s work better thus motivating the higher 
fraction of success of the Burkert profile. 

The bias values TZ(x) are summarized in Table 5 for the 
case with 5=1/3 and shows that, even assuming a Burk- 
ert halo, the disc mass turns out to be overestimated hence 
mimicking maximal disc solutions. On the contrary, the halo 
virial mass is only modestly biased, especially if compared 
to the outcome of previous fits. Although the lZ(M v i r ) dis- 
tribution has a long tail to the right of the mean value, 
H(M V i r ) ^ 2 for most of the galaxies in the WS sample. 
Since Jdm depends on 1/Md and Md is overestimated (while 
the dependence on Mh(Ropt) is weaker being this quantity 
present both at the numerator and denominator), it is then 
expected that Jdm is biased low. This is indeed what we 
find in accordance with the similar result obtained for the 
PI model. Note, however, that this time the bias is smaller 
than for the PI model, even if still quite significant. 

Finally, we have repeated the above analysis setting 8 = 
1.0 when simulating the rotation curves. The WS sample 
still contains 73% of the full simulated sample, but only 
17% of them pass the criteria for entering the BS sample. 
This situation is quite similar to what takes place for the PI 
model so that we do not discuss it here anymore. 



6 CONCLUSIONS 

General Relativity has been experimentally tested on scales 
up to the Solar System one so that assuming it still holds 
on much larger scales, such as the galactic and cosmological 
ones, is actually nothing else but an extrapolation. Moti- 
vated by this consideration and the difficulties in explaining 
the observed accelerated expansion without introducing new 
unknown ingredients like dark energy, a great interest has 
been recently devoted to modified gravity theories which 
have proven to work remarkably well in fitting the data and 
predicting the correct growth of structures. Modifying Gen- 
eral Relativity has impact at all scales so that, provided no 
departures from standard results, well established at Solar 
System scales, one cannot exclude a priori that the gravita- 
tional potential generated by a point mass source has not 
the usual Keplerian fall off, tf> oc 1/r, but a weaker one. Here 
we have considered the case of a Yukawa- like correction, 
i.e. <f> oc (l/r)[l + <5exp(— r/A)] where the scale length A is 
related to the effective scalar field (coupled with matter) in- 
troduced by several modified theory of gravity. In particular, 
this kind of potential comes out in the weak field limit of 
/(i?)-gravity. Provided A is much larger than the Solar Sys- 
tem scale, the corrections to the potential can significantly 
boost the circular velocity for an extended system like a 
spiral galaxy. Assuming Newtonian gravity to compute the 
theoretical rotation curve to the observed one then intro- 
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duces a systematic error which can bias the estimate of the 
galaxy parameters thus leading to misleading conclusions. 

To investigate this issue, we have fitted the Newtonian 
circular velocity of some widely used halo models to a large 
sample of simulated rotation curves computed using the as- 
sumed modified potential and a disc + NFW profile. Com- 
paring the input parameters with the best fit ones allows us 
to draw some interesting lessons on the consequences of a 
systematic error in the adopted gravity theory. As a general 
result, we find that, for cusped halo models, the disc mass is 
underestimated, while the halo scalelength and virial mass 
are biased high. Since the concentration c is also biased low, 
the c- M v i r relation will have a different slope and intercept 
than the one we have used to generate the model based on 
the outcome of N- body simulations. Moreover, if left free in 
the fit, the inner slope of the density profile turns out to be 
biased low if A is smaller than the disc half mass radius Rd- 
It is interesting to note that halo models shallower than the 
NFW one in the inner regions with (c, M v i r ) values not con- 
sistent with the prediction of N - body simulations and with 
maximal discs are a common outcome of fitting observed 
(not simulated) rotation curves. Such results are usually in- 
terpreted as failures of the ACDM model due to misleading 
assumptions on the dark matter particles properties (such 
as their being hot or cold and their interaction cross sec- 
tions) or to having neglected the physics of baryons in the 
simulations. Our analysis point towards a different expla- 
nation considering these inconsistencies as the outcome of 
forcing the gravitational potential to be Newtonian when 
it is not. In order to test such an hypothesis, one should 
fit a homogenous set of well sampled and radially extended 
rotation curves assuming a theoretically motivated modi- 
fied gravity potential (to set the strength S of the corrective 
term) and a classical NFW model (since it is in agreement 
with N-body simulations also in fourth order theories). 

Should the scalelength A of the modified potential be 
smaller than the disc one Rd, we have shown that forcing 
the potential to be Newtonian may also lead to completely 
mismatch not only the model parameters, but also the halo 
density profiles. Indeed, cored models, such as the PI and 
Burkert ones, turn out to be able to fit equally well the 
rotation curve data. We have not carried out here a case- 
by-case comparison among the different models considered 
since this will depend critically on what the actual (not the 
simulated) uncertainties are and on the sampling and ra- 
dial extent of the data. We can however anticipate that, 
for most cases, the cored models will work better than the 
cusped ones since they are better able to increase the in- 
ner rotation curve redistributing the total dark matter mass 
inside the core thus leaving almost unchanged the outer ro- 
tation curve. Such a result may have deep implications on 
the cusp/core controversy which should then be read as an 
evidence of an inconsistent assumption about the correct un- 
derlying gravity theory. From an observational point of view, 
one could try to fit our NFW + modified potential to LSB 
galaxies since this dark matter dominated systems are usu- 
ally considered the best examples of the cusp/core problem. 
Note that our simulated sample does not actually contain 
LSB - like galaxies since we have set the input dark matter 
mass fraction as /dm ~ 50% and adjusted the disc parame- 
ters having in mind a Milky Way - like spiral galaxy. Should 
we have simulated only LSB - like systems, we expect that 



cored models will definitely be preferred over cusped ones 
thus further strengthening our conclusions. 

As a general remark, we have also obtained that it is 
actually quite difficult (if not impossible) to fit in a satisfac- 
tory way the simulated rotation curves with Newtonian halo 
models if the amplitude of the modified potential is set to 
5=1, i.e. when the Yukawa term in the point mass case has 
the same weight as the Keplerian one. This is not surprising 
since the boost in the circular velocity in the halo dominated 
regions is so large that can only be reproduced by increas- 
ing the virial mass to unacceptably large values. Although a 
more detailed analysis is needed, this result could be consid- 
ered as an evidence against a too large deviation from the 
Keplerian 1/r scaling. Indeed, should the case 5 = 1.0 be 
realistic, then the actually observed rotation curves would 
resemble the simulated ones and hence we should have been 
unable to fit them. This is obviously not the case since a 
plethora of successful fits are available in literature. We can 
therefore argue that the case S — 1 is not realistic at all or, 
in other words, that a Yukawa - like deviation from the New- 
tonian potential must be only a subdominant correction to 
the 1/r scaling in the inner galaxy regions. 

Although a more detailed analysis is needed, we would 
finally stress that our analysis points towards a new usage of 
the rotation curves. Looking for inconsistencies rather than 
for agreement between these data and Newtonian models 
can indeed tell us not only whether the dark matter particles 
properties should be modified or not, but also whether our 
assumptions on the underlying theory of gravity are correct 
or not. Although it is likely that a definitive answer on this 
question could not be achieved in this way, the analysis of the 
rotation curves data stands out as a new tool to deal with 
modified gravity at scales complementary to those tested 
by cosmological probes. Asking for consistency among the 
results on such different scales could help us to select the 
correct law governing the dominant force of the universe. 
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